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1.  Introduction: 


The  broad  objective  of  the  proposed  research  is  to  develop  an  optimized  system  design  and 
associated  image  reconstruction  algorithms  for  a  hybrid  three-dimensional  (3D)  breast  imaging 
system  that  combines  OAT  and  UST.  I  have  made  excellent  progress  during  last  year  in 
accomplishing  the  specific  tasks  for  this  project.  In  the  past  year,  my  research  has  primarily  been 
focused  on  (i)  developing  time-of-flight  extraction  algorithms  to  perform  USCT,  (ii)  developing 
image  reconstruction  algorithms  for  USCT,  (iii)  developing  OAT  algorithms  (iv)  accelerating 
OAT  algorithm  to  enable  3D  image  reconstruction  for  breast  imaging,  (v)  evaluating  and 
validating  algorithms  by  computer  simulation  studies  and  experimental  phantom  studies. 

2.  Keywords: 


Ultrasound  imaging,  optoacoustic  imaging,  photoacoustic  imaging,  iterative  image 
reconstruction,  nonlinear  optimization,  breast  imaging,  GPU  acceleration 

3.  Overall  Project  Summary 

Task  1:  Construct  a  computational  model  of  the  OAT/UST  imager  and  identify  optimal 
system  geometries: 

Computational  modeling  of  the  proposed  imager: 

Imager  development  is  based  on  comprehensive  computer  models  of  OAT  and  UST.  The  system 
design  studies  are  conducted  concurrently  with  the  development  of  the  image  reconstruction 
algorithms  so  that  they  can  be  informed  and  refined  jointly. 

Optimization  studies: 

I  have  conducted  numerical  studies  to  obtain  an  optimal  imager  design.  The  ultrasound 
tomographic  image  quality  depends  strongly  on  the  distribution  and  number  of  emitters  and 
transducers  pairs.  Moreover,  image  quality  also  depends  on  the  number  of  tomographic  views. 
These  parameters  have  to  be  carefully  chosen  to  obtain  a  feasible  experimental  design. 
Optimization  studies  were  performed  to  detennine  the  number  of  emitters  and  tomographic 
views  for  a  fixed  receiver’s  array.  To  obtain  3D  SOS  distribution,  the  emitters  will  be 
distribution  on  a  planar  surface. 


Task  2:  Development  of  reconstruction  methods  for  sparse-array  3D  UST 
Reconstruction  of  SOS  distribution: 

Iterative  image  reconstruction  algorithms  have  been  developed  to  reconstruct  the  SOS 
distribution.  These  algorithms  have  also  been  validated  for  experimental  phantom  studies.  To 
perform  USCT,  measured  time-of-flight  (TOF)  was  extracted  from  the  measured  signals  of  the 
transducer  elements. 

Time-of-flight  (TOF)  extraction  of  the  transmission  ultrasound  signals:  I  utilized 
geometrical  acoustic-based  ray  theory  to  establish  a  non-linear  model  that  relates  the  measured 
TOF  values  to  speed  of  sound  (SOS)  distribution.  To  solve  this  nonlinear  optimization  problem, 
we  needed  to  extract  time-of-flight  (TOF)  from  the  measured  signal.  This  is  a  very  important 
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pre-processing  step  for  good  image  quality. 
In  search  of  the  best  TOF  extraction 
technique,  six  TOF  extraction  algorithms 
have  been  implemented  and  compared:  (i) 
envelope-detection  method,  (ii)  picking  the 
max  value  of  filtered  signal,  (iii)  AIC- 
method,  (iv)  weighted-AIC  picker,  (v)  cross¬ 
correlation  method,  and  (vi)  thresholding 
method  from  windowed  and  filtered  signal 
[1].  These  methods  were  investigated  for 
both  the  computer  simulation  by  adding 
different  levels  of  noise  and  the  experimental 
data  (provided  by  Tomowave  Laboratories 
Inc.,  Houston  TX). 

Bent-Ray  method:  Algorithms  are 

developed  for  reconstructing  the  SOS 
distribution  of  breast  from  knowledge  of 
time-of-flight  (TOF)  measurements  of  the 
transmission  ultrasound  signals.  Utilizing  the 
geometrical  acoustic-based  ray  theory,  a  non- 


-  Phantom 

-  flay  method 


y  slice 


Figure  1:  (a)  Numerical  breast  phantom;  (b)  image 
reconstructed  using  adjoint-state  method;  (c)  image 
reconstructed  using  bent-ray  algorithm;  (d)  plot 
corresponding  to  vertical  line  in  (a);  (e)  plot  corresponding 
to  horizontal  line  in  (a).  To  have  a  fair  comparison  between 
two  methods  no  regularization  has  been  used  in  both 
methods. 


linear  model  has  been  established  that  relates  the  measured  TOF  values  to  the  SOS  distribution. 


For  a  given  SOS  distribution,  numerically  solution  of  the  Eikonal  equation  yielded  the  ray  paths. 
An  iterative  reconstruction  method  was  developed  for  inverting  the  resulting  system  of  equations 
that  alternatively  updates  the  estimates  of  the  SOS  and  ray  paths,  minimizes  a  regularized  cost 
function  to  obtain  the  final  estimate  of  the  SOS  [2]. 


Adjoint-State  Method:  I  also  investigated  a 
partial  differential  equation-based  Eulerian 
approach  to  travel-time  tomography  as  an 
alternative  approach  [3].  The  work  on 
comparison  of  the  Adjoint-State  Method  are 
Ray-tracing  algorithm  was  presented  in  SPIE 
Photonics  West,  2014  and  SPIE  Medical 
Imaging  2014.  For  detail  implementation  of 
this  algorithm  please  see  attached  proceedings 
paper.  The  wavefonn  inversion  method  for 
SOS  reconstruction  has  also  been  explored  in 
the  group.  In  this  regard,  Adjoint-State  method 
provides  a  suitable  initial  SOS  distribution  to 
aid  waveform  inversion  method. 


Tumor 

Diameter 

(mm) 

X 

(FWHM) 

(mm) 

Y 

(FWHM) 

(mm) 

SOS 

(Actual 

1600 

m/s) 

Error 

(%) 

12 

12 

11.75 

1598.2 

0.11 

10 

9.75 

9.5 

1598.3 

0.1063 

8 

7.75 

7.25 

1603.1 

-0.194 

6 

7.0 

6.75 

1582.9 

1.068 

4 

6.5 

6 

1560.0 

2.5 

Table  1:  Performance  of  the  adjoint-state  method  has  been 
evaluated  be  performing  many  numerical  simulations 
studies  very  varying  tumor  sizes.  The  adjoint-state  method  - 
--  performs  well  for  smooth  SOS  distributions  —  can  give 
accurate  estimate  of  the  size  and  SOS  value  for  tumors 
ranging  in  6  mm  to  12  mm  tumor  sizes. 


I  perfonned  several  numerical  studies  to 
compare  bent-ray  and  adjoint-state  method. 

Figure  (1)  shows  the  comparison  of  the  two  methods  for  a  numerical  breast  phantom.  This 
numerical  study  was  performed  for  a  ring  scanner  consisting  of  256  transducer  elements.  In 
another  study,  I  varied  the  tumor  size  and  performed  image  reconstruction  to  access  tumor 
detectability  using  adjoint-state  method.  Table  (1)  summarizes  results  from  this  study.  It  can  be 
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Figure  2:  (a)  Delta-TOF  for  the  experimental  three-tube  phantom;  (b) 
USCT  image  reconstructed  using  adjoint-state  method;  (c)  OAT  image 


seen  that  the  tumors  with  diameter  large  than  4  mm  can  be  characterized  using  adjoint-state 
method. 

The  method  was  also  validated  for  the  experimental  phantom  studies.  In  this  study,  the  pressure 
signals  were  recorded  using  a  transducer  array  consisting  of  64  elements  and  a  single  Laser 
Induced  Source  placed  opposite  to  the  middle  transducer.  The  data  was  recoded  for  150  views. 
The  experimental  phantom  consisted  of  three  tubes  at  varying  salt  concentration  to  exhibit 
different  acoustic  and  optical  properties.  Results  are  shown  in  Figure  (2)  for  both  the  SOS 
distribution  and  optical  absorbed  energy  density. 


Reconstruction  of  reflectivity:  An  algorithm  has  been  developed  to  produce  reflectivity  maps. 
These  algorithms  are  based  on  the  Synthetic  Transmit  Aperture  (STA)  approach  [4].  This  method 
utilizes  multiple  elements  or  a  single  source  to  produce  the  spherical  waves  and  the  whole  image 
is  being  reconstructed  for  each  emitted  signal.  The  final  reflectivity  map  is  obtained  by 
accumulating  these  individual  images.  It  has  been  shown  that  the  STA  method  improves  SNR. 
This  method  is  especially  useful  in  our  current  study  because  we  are  using  laser  induced 
ultrasound  emitter  (LUS),  which  produce  spherical  waves. 


Task  3:  Ultrasound-assisted  OAT  image  reconstruction: 

Development  of  imaging  models  and  reconstruction  algorithms:  An  interpolation-based 
discrete-discrete  imaging  model  has  been  implemented  to  perfonn  3D  OAT  for  breast  imaging 
[5].  In  the  new  implementation,  an  unmatched  back  projection  (or  pixel-driven)  scheme  has  been 
used  and  validated  in  computer  simulations  studies.  This  algorithm  is  five  times  more  efficient 
than  the  ray-driven  back-projection  and  allows  to  perfonn  iterative  image  reconstruction  for 
large  fields-of-views,  making  it  very  suitable  for  breast  imaging.  To  efficiently  mitigate  data 
incompleteness,  noise,  and  model  error,  I  investigated  the  least-squares  objective  regularized  by 
a  TV-nonn  penalty.  I  implemented  the  fast  iterative  shrinkage/thresholding  algorithm  (FISTA)  to 
minimize  cost  function  with  TV  regularization  [6]. 

GPU  implementation  of  image  reconstruction  algorithms:  Improved  GPU-based 
implementations  of  a  numerical  imaging  model  and  its  adjoint  have  been  developed  for  use  with 
general  gradient-based  iterative  image  reconstruction  algorithms.  Particularly,  two  types  of 
computation-reduced  discretization  methods  have  been  employed;  a  parallel  fast  GPU-based 
Fourier  transform  (FFT)  algorithm  was  employed  to  accelerate  the  calculation  of  the  temporal 
convolution  with  ultrasonic  transducer  responses;  and  a  volume-reduction  method  is  proposed  to 
reduce  the  computation  for  applications  with  irregular  field-of-view  (breast  imaging).  The  results 
suggest  that  the  proposed  implementation  is  more  than  five  times  faster  than  previous 
implementations  for  a  single  GPU.  In  addition,  the  algorithm  has  also  been  developed  to  use 
multiple  GPUs  further  reducing  the  computational  time.  The  work  will  be  presented  in  SPIE 
Photonics  West,  2015. 
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Task  4:  Validate  prototype  imager  and  image  reconstruction  algorithms 


Phantom  imaging  studies:  The  imager  and 
algorithm  designs  have  been  infonned  and 
evaluated  by  use  of  experimental  studies  for 
well-characterized  multi-modality  phantoms 
conducted  at  TomoWave  Laboratories  under 
Dr.  Oraevsky.  We  used  phantoms  that  have 
tumors  located  at  different  depths  and  have 
different  optical  absorption  properties  to 
quantify  the  sensitivity  of  the  OAT  system. 
Simplified  versions  of  the  phantoms  will  be 
imaged  for  characterizing  the  spatially 
variant  spatial  resolution  and  noise  properties 
of  the  reconstructed  images. 


Figure  3:  (a)  LOUIS-3DB  imaging  module,  (b)  3D  rendered 
reconstructed  image  from  experimental  phantom  study.  The 
image  is  obtained  using  accelerated  iterative  image 
reconstruction  algorithm.  Both  tumor  inserts  and  blood  vessel 
phantom  are  clearly  visible  in  the  reconstructed  image. 


4.  Key  Research  Accomplishments: 

Accelerating  three-dimensional  iterative  image  reconstruction  algorithm:  One  of  the  key 
research  accomplishment  was  the  successful  development  of  accelerated  iterative  image 
reconstruction  algorithm  for  3D  OAT  imaging.  The  task  to  perform  3D  OAT  breast  imaging 
presents  many  challenges.  One 
of  the  challenge  is  to  account 
for  large  field-of-view.  In  a 
typical  study,  image 
reconstruction  volume  is  120  x 
120  x  90  mm3  and  the  existing 
iterative  image-reconstruction 
algorithms  present  non-practical 
image  reconstruction  durations 
(days  to  obtain  a  single  high- 
resolution  image).  Analytical 
algorithms,  are  efficient  e.g. 
filter-back  project  (FBP),  cannot 
account  for  data  inconsistencies, 
model  error,  and  noise  in  the 
data.  This  makes  the  accelerated 

image  reconstruction  algorithm  very  important  milestone  of  the  proposed  project  and  will  also 
benefit  the  optoacoustic  tomography  research  field.  Figure  (4)  shows  results  for  a  numerical 
phantom  study  for  the  accelerated  interpolation-based  imaging  model. 


Figure  4:  Slices  for  (a)  original  phantom  (b)  reconstructed  image  using  the 
original  interpolation-based  algorithm  with  matched  back  projection 
scheme;  (c)  reconstruction  image  using  accelerated  algorithm,  which  is  20- 
times  faster  and  exhibits  same  image  quality;  (d)  line  plot  corresponding  to 
vertical  line  in  (a);  (e)  profile  for  horizontal  line  in  (a)  for  (a)  (solid  line),  (b) 
(dashed  line)  and  (c)  (dotted  line)  respectively. 


5.  Conclusion 

I  will  continue  to  improve  efficiency  and  accuracy  of  the  reconstruction  algorithms  for  USCT 
and  OAT  to  perfonn  breast  imaging.  Algorithms  to  reconstruct  attenuation  will  also  be 
developed.  I  will  perform  image  reconstruction  for  clinical  studies  as  soon  as  clinical  data 
becomes  available  and  will  use  the  experimental  study  to  investigate  and  improve  the  imaging 
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algorithms  and  imager  design.  Future  studies  will  also  be  focused  on  developing  task-based 
optimization  studies  for  OAT  and  USCT  next  generation  imager  design.  I  will  also  develop 
algorithms  to  perform  USCT-assisted  OAT  imaging. 

6.  Publications,  Abstracts,  and  Presentations 

la.  Waveform  Inversion  with  Source  Encoding  for  Breast  Speed-of-Sound  Reconstruction  in 
Ultrasound  Computed  Tomography  (submitted) 

Authors:  Kun  Wang,  Thomas  Mathew,  Fatima  Anis,  Cuiping  Li,  Neb  Duric,  and  Mark  Anastasio 
Journal:  IEEE  Transaction  on  Medical  Imaging 

lb*.  Investigation  of  the  adjoint-state  method  for  ultrasound  computed  tomography:  a  numerical 
and  experimental  study 

Authors:  Fatima  Anis,  Yang  Lou,  Andre  Conjusteau,  Sergey  Ermilov,  Alexander  Oraevsky  and 
Mark  A.  Anastasio 

Conference:  SPIE  Photonics  West-  2014,  San  Francisco  CA 

2b-  Investigation  of  a  method  for  laser-induced  ultrasound  tomography  that  eliminates  the  need 
for  ray-tracing 

Authors:  Fatima  Anis,  Yang  Lou,  Andre  Conjusteau,  Sergey  Ermilov, 

Alexander  Oraevsky  and  Mark  A.  Anastasio 
Conference:  SPIE  Medical  Physics-  2014,  San  Diego  CA 

3b-  Title:  Accelerated  iterative  image  reconstruction  in  three-dimensional  optoacoustic 
tomography 

Authors:  Fatima  Anis,  Yang  Lou,  Kun  Wang,  Richard  Su,  Tanmayi  Oruganti,  Andre 
Conjusteau,  Sergey  Ermilov,  Alexander  A.  Oraevsky  and  Mark  A.  Anastasio 
Conference:  SPIE  Photonics  West,  2015 

4b-  Title:  Waveform  Inversion  with  Source  Encoding  for  Breast  Speed-of-Sound  Reconstruction 
in  Ultrasound  Computed  Tomography 

Authors:  Kun  Wang,  Thomas  Mathew,  Fatima  Anis,  Cuiping  Li,  Neb  Duric,  and  Mark 
Anastasio 

Conference:  SPIE  Medical  Physics,  2015 

7.  Training 

The  year  has  been  very  fruitful  as  I  continue  to  benefit  from  many  scholarly  activities  and  kept 
adding  to  my  skills  as  imaging  scientist. 

1-  Attending  course  (E62  BME  500  67):  In  Fall  2014, 1  participated  in  BME  course  on  the 
imaging  science.  The  class  held  weekly  for  two  hours.  Some  of  the  study  topics  included 
continuous  and  discrete  object  representations,  imaging  operators,  image  statistics, 
imaging  quality  assessment,  ideal  observer,  Hotelling  observers  and  imaging  errors.  The 
course  provided  me  with  a  formal  training  as  biomedical  image  scientist. 

2-  Conferences: 

i-  SPIE  Photonics  West,  2014  San  Francisco  CA:  The  conference  is  the  major 
international  meeting  held  annually  for  the  optoacoustic/photoacpustic  imaging.  I 
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attended  this  meeting  and  was  greatly  benefited  from  the  presentation  and  poster  sessions 
as  well  as  constructive  meetings  with  the  Prof.  Oraevsky  and  his  team  about  the  hybrid 
OAT/  USCT  breast  imager. 

ii-  SPIE  Medical  Imaging,  2014,  San  Diago  CA:  This  conference  is  another  major 
international  conference,  which  holds  annually  and  covers  broad  range  of  topics 
concerning  medical  imaging  and  diagnostics.  I  was  specially  benefited  from  many  talks 
and  poster  presentations  about  the  image  quality  assessment.  Moreover,  the  dedicated 
Ultrasonic  Imaging  and  Tomography  sessions  on  ultrasound  imaging  provided  me  with  a 
great  opportunity  to  learn  about  ultrasound  tomography  in  medical  imaging. 

3-  Visiting  TomoWave  Inc,  Houston  TX: 

I  visited  TomoWave  Inc.  in  December  2013.  During  the  visit,  I  benefited  from  lab  tours 
and  learned  about  practical  aspects  of  OAT/USCT  imaging. 

4-  Algorithm  development  and  GPU  computing:  I  continued  to  establish  more  skills 
towards  algorithm  development.  One  of  the  major  achievements  towards  this  end  was  to 
leam  CUDA  programming  from  other  group  members.  This  training  will  continue  to 
benefit  me  through  the  remainder  of  the  project. 

8.  Reportable  Outcomes 

Nothing  to  report 

9.  Other  Achievements 

Nothing  to  report 
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ABSTRACT 

In  this  work,  we  investigate  a  novel  reconstruction  method  for  laser-induced 
ultrasound  tomography  (UST)  breast  imaging  that  circumvents  limitations  of 
existing  methods  that  rely  on  ray-tracing.  There  is  currently  great  interest  in 
developing  hybrid  imaging  systems  that  combine  optoacoustic  tomography 
(OAT)  and  UST.  There  are  two  primary  motivations  for  this:  (1)  the  speed- 
of-sound  (SOS)  distribution  reconstructed  by  UST  can  provide 
complementary  diagnostic  information;  and  (2)  the  reconstructed  SOS 
distribution  can  be  incorporated  in  the  OAT  reconstruction  algorithm  to 
improve  OAT  image  quality.  However,  image  reconstruction  in  UST 
remains  challenging.  The  majority  of  existing  approaches  for  UST  breast 
imaging  involve  ray-tracing  to  establish  the  imaging  operator.  This  process 
is  cumbersome  and  can  lead  to  severe  inaccuracies  in  the  reconstructed  SOS 
images  in  the  presence  of  multiple  ray-paths  and/or  shadow  zones. 

To  circumvent  these  problems,  we  implemented  a  partial  differential 
equation-based  Eulerian  approach  to  UST  that  was  proposed  in  the 
mathematics  literature  but  never  investigated  for  medical  imaging 
applications.  This  method  operates  by  directly  inverting  the  Eikonal 
equation  without  ray-tracing.  A  numerical  implementation  of  this  method 
was  developed  and  systematically  compared  to  existing  reconstruction 
methods  for  UST  breast  imaging.  We  demonstrated  the  ability  of  the  new 
method  to  reconstruct  accurate  SOS  maps  from  TOF  data  obtained  by  a  3D 
hybrid  OAT/UST  imager  built  by  our  team. 
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ABSTRACT 

In  this  work,  we  investigate  a  novel  reconstruction  method  for  laser-induced  ultrasound  computed  tomography 
(USCT)  breast  imaging  that  circumvents  limitations  of  existing  methods  that  rely  on  ray-tracing.  There  is 
currently  great  interest  in  developing  hybrid  imaging  systems  that  combine  optoacoustic  tomography  (OAT)  and 
USCT.  There  are  two  primary  motivations  for  this:  (1)  the  speed-of-sound  (SOS)  distribution  reconstructed  by 
USCT  can  provide  complementary  diagnostic  information;  and  (2)  the  reconstructed  SOS  distribution  can  be 
incorporated  in  the  OAT  reconstruction  algorithm  to  improve  OAT  image  quality.  However,  image  reconstruction 
in  USCT  remains  challenging.  The  majority  of  existing  approaches  for  USCT  breast  imaging  involve  ray¬ 
tracing  to  establish  the  imaging  operator.  This  process  is  cumbersome  and  can  lead  to  inaccuracies  in  the 
reconstructed  SOS  images  in  the  presence  of  multiple  ray-paths  and/or  shadow  zones.  To  circumvent  these 
problems,  we  implemented  a  partial  differential  equation-based  Eulerian  approach  to  USCT  that  was  proposed 
in  the  mathematics  literature  but  never  investigated  for  medical  imaging  applications.  This  method  operates 
by  directly  inverting  the  Eikonal  equation  without  ray-tracing.  A  numerical  implementation  of  this  method  was 
developed  and  compared  to  existing  reconstruction  methods  for  USCT  breast  imaging.  We  demonstrated  the 
ability  of  the  new  method  to  reconstruct  SOS  maps  from  TOF  data  obtained  by  a  hybrid  OAT/USCT  imager 
built  by  our  team. 

Keywords:  ultrasound  tomography,  optoacoustic  tomography,  photoacoustic  tomography,  breast  cancer  imag¬ 
ing 


1.  INTRODUCTION 

Transmission  ultrasound  computed  tomography  (USCT)  is  an  emerging  imaging  modality  with  many  biomedical 
applications.  USCT  can  be  employed  to  retrieve  anatomical  information  of  tissues  e.  g.  speed  of  sound,  acoustical 
impedance  and  reflectivity.  The  effectiveness  of  USCT  in  tumor  detection  has  been  discussed  in  recent  studies. 1-3 
It  is  known  that  cancerous  tissues  have  higher  SOS  values  compared  to  the  benign  fatty  masses  and  healthy 
breast  tissues.  A  clinical  ultrasound  ring  array  scanner  for  breast  cancer  diagnosis  (Computed  Ultrasound  Risk 
Evaluation  (CURE))  has  been  proposed.1,2  This  system  consists  of  256  transducer  elements  distributed  on 
a  ring  with  a  20  cm  diameter.  Another  prototype  has  also  been  developed  by  SoftVue  and  consists  of  2048 
transducers.4  The  system  is  capable  of  reconstructing  a  series  of  2D  slices  of  the  SOS,  acoustic  attenuation,  and 
reflectivity  distributions.  Techniscan  (Salt  Lake  City  UT)  introduced  a  commercial  USCT  system  that  employs 
three  transducer  probes  placed  around  the  breast.  The  transducer  system  is  mechanically  rotated  to  reconstruct 
2D  slices  and  subsequently  vertically  scanned  to  capture  multiple  slices  to  obtain  3D  images.5  Finally,  an 
ultrasound  imaging  module  capable  of  generating  three-dimensional  SOS  distributions  has  been  investigated  by 
researchers  at  the  Karlsruhe  Institute  of  Technology  (KIT),  Germany.6 

Biomedical  applications  of  USCT  commonly  employ  geometrical  acoustics  models  and  require  time-of-flight 
(TOF)  measurements.  TOF  data  that  have  been  recorded  for  many  source-receiver  pairs  can  be  employed  for 
reconstruction  of  the  SOS  distribution.  The  reconstruction  of  the  speed  of  sound  distribution  is  conventionally 
performed  by  using  ray-tracing  (RT)  methods.1,2,7,8  To  account  for  the  curvature  of  the  ray  paths,  the  rays  are 
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traced  along  the  negative  gradient  of  the  TOF  distribution.8  Ray-tracing  can  become  cumbersome,  especially 
for  three-dimensional  USCT.  Moreover,  unlike  X-ray  computed  tomography,  the  heterogeneous  SOS  distribution 
results  in  uneven  ray-distributions,  which  makes  the  inverse  problem  ill-conditioned.  In  this  work,  we  will 
investigate  a  different  approach  for  SOS  image  reconstruction.  This  method,  the  adjoint  state  (AS)  method,  has 
previously  been  employed  for  seismic  tomography.9, 10  We  are  investigating  the  AS  method  for  USCT  for  the 
first  time  for  biomedical  applications. 


2.  DESCRIPTION  OF  NUMERICAL  STUDIES 

To  perform  USCT  using  RT,  we  have  employed  the  the  geometrical  ray  theory  for  sound  waves.  This  approxi¬ 
mation  results  in  a  non-linear  model,  the  eikonal  equation,  to  relate  the  measured  TOF  values  to  the  speed  of 
sound  distribution  as 

|VT(r)|  =  -L.  (1) 

c(r) 

In  Eq.  (1),  VT  is  the  gradient  of  the  TOF,  T,  and  c  is  the  SOS  distribution,  both  of  which  are  a  function  of 
position  as  denoted  by  the  position  vector  rdl2.  Currently,  the  bent-ray  reconstruction  is  a  widely  employed 
reconstruction  technique  for  USCT  because  it  incorporates  refraction  during  sound  wave  propagation.  The 
eikonal  equation  is  solved  numerically  by  finite  difference  methods11  to  obtain  a  TOF  map  for  a  given  source 
corresponding  to  a  certain  speed  of  sound  map  c(r). 


2.1  Ray-tracing  reconstruction  method 

In  RT  methods,  the  TOF  is  calculated  as  the  line-integral  over  the  slowness  distribution  over  the  ray-path 
connecting  the  source  and  the  receiver  location: 


T(r) 


(2) 


The  dependence  of  the  ray-path,  T(c),  on  the  SOS  distribution  makes  it  a  non-linear  problem.  The  discretized 
imaging  model  is  given  by 


T  =  H(c)  — ,  (3) 

c 

where  T  is  a  vector  of  TOF  measurments,  c  is  a  finite-dimensional  representation  of  the  SOS,  and  H(c)  is  the 
system  matrix.  To  formulate  H(c),  we  implemented  a  RT  method.  Weights  were  assigned  to  pixels  in  the  discrete 
SOS  map  based  on  the  number  of  times  each  pixel  was  intersected  by  the  rays.  This  weight  matrix  constitutes 
H(c). 

To  estimate  the  SOS  distribution  from  the  measured  TOF  data,  we  solved  the  following  optimization  problem: 

c  =  argmin  ||  T  —  T*  ||2  +vg(c),  (4) 

C 

where  c  denotes  the  sought-after  estimate  of  the  the  SOS  distribution,  T*  is  the  measured  TOF  data  from 
all  source-transducer  pairs,  g{ c)  is  a  penalty  function,  v  is  a  regularization  parameter,  and  T  is  the  computed 
TOF  found  by  solving  (1).  To  minimize  Eq.  (4)  we  used  the  Limited  BFGS  method.12  In  solving  the  nonlinear 
optimization  problem,  we  evaluated  the  gradient  of  Eq.  (4)  as 

Vc  =  2H(c)t  [H(c)c  —  T*]  +  vVg(c).  (5) 

It  should  be  noted  from  the  above  equation  that  the  first  term  is  a  linear  approximation  of  the  true  non-linear 
gradient  of  the  objective  function.  The  above  linearized  gradient  is  primarily  used  in  USCT  for  biomedical 
applications. 
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2.2  Adjoint-state-based  reconstruction  method 

We  implemented  a  previously  proposed  algorithm  for  USCT  reconstructed  based  on  the  adjoint-state  method.9 
The  mismatch  energy  functional  between  the  measured  and  simulated  data  is  defined  as9 

E[c(r)}  =  l  f  \T(r)-T*(r)fdn,  (6) 

J  s 

where  T*|s  is  the  measured  TOF  and  Tjg  is  computed  by  solving  Eq.  (1).  The  quantity  in  Eq.  (6)  —  the 
energy  functional  —  measures  the  L2-difference  between  the  the  solution  of  the  eikonal  equation,  T,  and  the 
experimental  measurement,  T*,  on  the  measurement  surface  S.  Using  the  adjoint-state  method  for  a  small 
perturbation  ec  to  c,  the  gradient  of  the  energy  is  defined: 


SE  =  e 


r)^(r) 

3(r) 


(7) 


Here,  V  is  volume  enclosed  the  measurement  surface  S  and  A(r)  is  the  adjoint  function  to  T(r)  that  satisfies  the 
following  adjoint  equation: 


V  •  [A(r)VT(r)]  =  0 


with  the  boundary  condition, 


[n.VT(r)]A(r)|s  =  [T*(r)-T(r)]s. 


(8) 

(9) 


Here  n  is  the  unit  outward  normal  of  the  surface  S.  To  minimize  the  energy  using  the  method  of  gradient 
descent,  a  perturbation  c(r)  =  —  A(r)/c3(r)  is  defined.  This  leads  to 


SE  =  — e  /  c2(r)dn  <  0,  (10) 

Jv 

where  V  denotes  the  region  interior  to  S.  By  solving  Eqs.  (8)  and  (9),  the  update,  c(r),  to  the  SOS  distribution 
can  be  computed.  Specifically,  the  SOS  distribution  is  updated  at  each  step  as: 


ck+ 1  =  ck  +  ekck 


(11) 


until  a  convergence  criterion  is  reached.  The  following  two  conditions  are  required  of  the  SOS  distribution: 
(i)  ck\s  =  0  and  (ii)  cfc+1  is  smooth.  To  fulfill  (ii),  a  regularization  term  similar  to  the  one  used  in  Eq.  (4) 
was  included.  The  filtering  scheme  defined  in  Lueng’s  work9  was  also  implemeted.  The  step  size,  ek,  can  be 
determined  by  using  the  Armijo-Golstein  rule  or  by  simply  setting  ek  =  e.  The  update  scheme  described  in 
Eq.  (11)  takes  a  large  number  of  iterations  to  converge.  Therefore,  we  used  the  limited-memory  Broydon- 
Fletclier  Goldfarb-Shanno  (L-BFGS)  method  to  solve  for  this  nonlinear  optimization  problem.  We  solved  the 
adjoint-state  equation  (8)  using  the  fast-sweeping  method9  with  the  boundary  conditions  defined  in  Eq.  (9). 


2.3  Experimental  Setup 

The  experimental  setup  consists  of  a  single  laser  ultrasound  (LU)  source  and  a  64  transducer  elements  arranged 
in  an  arc.  The  array  aperture  spans  a  152  degree  arc  with  a  radius  of  65  nun.  The  imaging  module  is  mounted 
and  centered  on  a  rotational  stage  operated  by  a  stepper  motor,  which  is  used  to  obtain  TOF  measurement  for 
150  views.  The  distance  between  the  central  element  of  the  arc  array  and  the  LU  source  is  130  mm.  Optoacoustic 
imaging  was  concurrently  performed.  More  detail  about  the  LU  sources  and  the  transducer  array  can  be  found 
elsewhere.13, 14 
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Figure  1.  (a)  Speed  of  sound  distribution  for  the  phantom;  reconstructed  SOS  distribution  for  the  RT  method  (b)  and  for 
the  AS  method  (c). 


3.  RESULTS  AND  DISCUSSION 

Computer-simulation  studies  were  conducted  to  compare  the  RT  method  with  the  AS  method.  In  this  study,  the 
objective  function  did  not  include  a  penalty  and  least  squares  estimates  of  the  SOS  distributions  were  computed. 
The  two-dimensional  phantom  depicting  the  SOS  distribution  is  shown  in  Figure  1(a).  The  phantom  consists 
of  three  discs  of  4.72  mm  diameter  with  constant  SOS  of  1.48,  1.6  and  1.7  mm//rs,  respectively.  The  pressure 
data  were  generated  using  the  k-Wave  software  package15  with  a  geometry  and  acoustic  properties  consistent 
with  our  experimental  system  design.  Gaussian  white  noise  was  added  to  the  calculated  pressure  signal  to  obtain 
experimentally-relevant  SNRs.  Figures  1(b)  and  1(c)  show  the  reconstructed  SOS  distribution  for  the  RT  method 
and  the  AS  method,  respectively.  In  the  case  of  the  RT  method,  streak  artifacts  are  very  visible  as  compared  to 
the  AS  method.  All  three  structures  can  be  seen  in  the  AS  reconstruction  of  the  SOS  distribution.  The  results 
show  that  the  AS  method  can  be  successfully  used  to  perform  USCT  for  biomedical  applications. 

To  check  the  accuracy  of  the  reconstructed  SOS  distribution,  we  selected  2  mm  x  2mm  regions  at  different 
locations  and  calculated  the  avaeraged  SOS  and  standard  deviation  in  those  region.  The  location  of  the  five 
regions  is  shown  in  Fig. 2(a).  The  bar  plot  in  Fig. 2(b)  shows  the  averaged  SOS  values  for  both  the  RT  and  AS 
method.  It  can  be  seen  that  AS  gives  accurate  SOS  values  for  the  selected  regions.  The  maximum  standard 
deviation  was  0.0226  nnn//is  for  region  ”A”. 

Finally,  we  studied  the  use  of  the  AS  method  to  perform  SOS  image  reconstruction  for  the  experimentally 
measured  TOF  from  our  LU  system.  The  experimental  phantom  consists  of  three  tubes  each  with  a  4.72  mm 
internal  diameter.  Tubes  were  filled  with  water  at  different  salt  concentrations  to  induce  different  SOS  values 
and  CUSO4  was  added  to  provide  optical  contrast.  For  this  case,  concurrent  optoacoustic  (OA)  and  ultrasonic 
data  acquisition  was  performed.  The  OA  reconstruction  of  the  phantom  is  shown  in  Fig.  3(a)  and  the  SOS 
reconstruction  is  shown  in  Fig.  3(b).  The  SOS  image  was  found  via  the  AS  reconstruction  method.  Once  again, 
all  three  discs  are  visible  in  the  reconstructed  SOS  distribution. 
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Figure  2.  (a)  Speed  of  sound  distribution  for  the  phantom  with  the  marked  location  of  five  regions;  (b)  bar  plot  for  the 
averaged  SOS  values  in  five  regions  for  phantom,  AS  method,  and  RT  method. 
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Figure  3.  (a)  Optoacoustic  image  of  three  tubes;  (b)  reconstructed  SOS  distribution  using  the  AS  method. 


4.  SUMMARY 

The  adjoint  state  method  has  been  implemented  for  biomedical  applications  of  USCT.  Images  reconstructed  from 
both  simulation  studies  and  measured  TOF  data  were  presented.  Ray  tracing  becomes  much  more  cumbersome 
for  three-dimensional  USCT.  Consequently,  the  adjoint  state  method  holds  great  promise  for  that  application. 
Further  numerical  studies  will  also  be  performed  to  quantify  resolution  and  noise  propagation  in  the  adjoint  state 
method. 
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Abstract 

Our  work  introduces  an  ultrasound  tomography  (UST)  reconstruction  algorithm  based  on  the  adjoint  method 
for  medical  imaging.  This  method  improves  current  ray-tracing  based  UST  reconstruction  algorithm  and  has 
been  previously  applied  to  seismic  travel-time  tomography  [S.  Leung  and  J.  Qian,  Comm.  Math.  Sci.  4 
(2006)].  Ultrasound  tomography  has  received  wide  attention  for  its  ability  to  help  breast  cancer  diagnosis 
both  by  providing  speed  of  sound  and  attenuation  information,  as  well  as  providing  adjunct  imaging  data  for 
optoacoustic  tomography  (OAT).  Current  image  reconstruction  algorithms  for  UST  are  usually  based  on  ray¬ 
tracing  and  gradient  methods.  Our  investigation  shows  two  drawbacks  of  these  methods  that  lead  to  inaccuracy 
in  image  reconstruction.  First,  ray  bending  in  ray-tracing  will  cause  an  uneven  distribution  of  ray  paths.  This 
will  lead  to  insufficient  updates  in  shadow  zones  (regions  covered  by  few  ray  paths)  and  cause  artifacts.  While 
this  effect  can  be  compensated  by  regularization  to  some  extent,  we  show  that  it  cannot  be  avoided  completely 
in  ray-tracing  methods.  Second,  often,  a  linear  approximation  of  the  gradient  objective  function  is  used,  which 
also  introduces  errors  into  the  gradient  descent  optimization  method.  We  will  demonstrate  that  using  the  adjoint 
method  to  directly  compute  the  Frechet  derivative  of  the  continuous  non-linear  objective  function  can  circumvent 
both  drawbacks  of  the  ray-tracing  method.  Numerical  simulations  are  then  given  to  show  the  improvement  of 
our  method  over  the  ray-tracing  method. 


1.  DESCRIPTION  OF  PURPOSE 

Transmission  ultrasound  tomography  (UST)  is  an  emerging  modality  that  has  a  spectrum  of  biomedical  ap¬ 
plications.  The  transmitted  ultrasound  signal  carries  anatomical  information  about  the  object,  e.  g.  speed  of 
sound,  acoustical  impedance  and  reflectivity.  Biomedical  applications  of  UST  commonly  employ  geometrical 
acoustics  models  and  require  time-of-ffight  (TOF)  measurements.  TOF  data  that  have  been  recorded  for  many 
source-receiver  pairs  can  be  employed  for  reconstruction  of  the  speed-of-sound  (SOS)  distribution. 

We  will  present  a  comparison  of  two  different  algorithms  for  the  accuracy  and  efficiency  to  perform  ultrasound 
computed  tomography.  One  of  the  methods  has  widely  been  utilized  for  the  medical  applications  of  ultrasound 
tomography.  In  this  ray-tracing  method  [Manohar,  et.  al.,  Appl.  Pliys.  Lett.  131911,  2007],  the  eikonal  equation 
is  solved  and  a  system  matrix  is  formulated  using  the  ray  paths  traced  from  the  receivers  locations  to  the  source 
locations.  The  other  method,  the  adjoint  -state  method  [S.  Leung  and  J.  Qian,  Comm.  Math.  Sci.  4  (2006)], 
eliminates  the  need  to  calculate  ray  paths  and  the  descent  direction  is  calculated  by  solving  the  adjoint  state 
equation. 
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2.  METHODS 


To  perform  UST,  we  have  used  the  the  geometrical  ray  theory  for  the  sound  waves.  This  approximation  results 
in  a  non-linear  model,  the  eikonal  equation,  to  relate  the  measured  time  of  flight  (TOF)  values  to  the  speed  of 
sound  distribution  (SOS) 


IWWI  =  ±  (1) 

In  Eq.  (1),  VT  is  the  gradient  of  the  travel  time,  T,  and  c  is  the  SOS  distribution,  which  both  are  a  function 
of  position  as  denoted  by  the  position  vector  r  G  R2.  Currently,  the  bent-ray  reconstruction  is  a  widely  used 
reconstruction  technique  for  UST  because  it  incorporates  refraction  during  sound  wave  propagation.  The  eikonal 
equation  is  solved  by  either  the  Fast  Marching  Method  or  a  Finite  Difference  Method  to  obtain  a  time-of-flight 
map  for  a  given  source  corresponding  to  a  certain  speed  of  sound  map  c(r). 

2.1  Ray-tracing-based  reconstruction  method 

We  write  the  optimization  problem  as: 


c  =  argmin  ||  T  —  T*  ||2  +vg(c),  (2) 

C 

where  c  is  the  discrete  representation  of  the  SOS  distribution,  c  is  the  estimate  of  the  the  SOS  distribution, 
T*  is  the  TOF  measurement  of  all  transducer  pairs,  g( c)  is  a  regularization  term  dependent  on  the  SOS,  v  is  a 
regularization  parameter,  and  T  is  the  computed  TOF  found  by  solving  (1).  Mathematically,  T  can  be  expressed 
as 

T  =  H(c)  — ,  (3) 

c 

where  H(c)  is  the  system  matrix.  To  formulate  H(c),  we  implemented  the  ray-tracing  method  between  each 
receiver  and  source  transducer  pair.  Weights  were  assigned  to  pixels  in  the  discrete  representation  of  the  object 
based  on  the  number  of  times  each  pixel  was  intersected  by  the  rays.  This  weight  matrix  constitutes  H(c). 
To  minimize  Eq.  (2)  we  used  the  Limited  BFGS  method.  In  solving  the  nonlinear  optimization  problem,  we 
evaluated  the  gradient  of  Eq.  (2)  as 


Vc  =  2H(c)t(H(c)c  —  T*)  +  v\/g(c).  (4) 

It  should  be  noted  from  the  above  equation  that  the  first  term  is  a  linear  approximation  of  the  true  non-linear 
gradient  of  the  objective  function.  The  above  linearized  gradient  is  primarily  used  in  UST  for  the  biomedical 
applications.  It  is  impractical  to  numerically  calculate  the  true  gradient  of  Eq.  (2). 

2.2  Adjoint-state-based  reconstruction  method 

We  define  the  mismatch  energy  functional  between  measured  and  simulated  data  as  [S.  Leung  and  J.  Qian, 
Comm.  Math.  Sci.  4  (2006)] 


E[c{v)]=l-[\T(v)-T*(v)\2dn,  (5) 

J  s 

where  c  is  the  speed  of  sound,  T*\$  is  the  measurement,  and  T\s  is  computed  by  solving  the  eikonal  equation, 
Eq.  (1).  The  quantity  in  Eq.  (5)  —  the  energy  functional  —  measures  the  L2-difference  between  the  the  solution 
of  the  eikonal  equation,  T,  and  the  experimental  measurement,  T* ,  on  the  measurement  surface  S.  Using  the 
adjoint-state  method  for  a  small  perturbation  ec  to  c,  we  define  the  gradient  of  the  energy: 


c(r)A(r) 

c3(r) 


(El 


SE  =  e 


(6) 


Here  A(r)  is  the  adjoint  variable  and  satisfies  the  following  adjoint  equation: 


V  •  [A(r)VT(r)]  =  0 


(7) 


with  the  boundary  condition, 


[n-VT(r)]A(r)|s  =  [r*(r)-T(r)]5.  (8) 

Here  n  is  the  unit  outward  normal  of  the  surface  S .  To  minimize  the  energy  using  the  method  of  gradient  descent 
we  choose  the  perturbation  c(r)  =  —  A(r)/c3(r).  This  leads  to 


5E  =  — e  f  c2(r)dn  <  0.  (9) 

Js 

By  solving  Eqs.  (7)  and  (8),  the  update  c(r)  to  the  SOS  distribution  can  be  computed.  Specifically,  the  SOS 
distribution  is  updated  at  each  step  by: 


ck+l  =ck  +  ek~k 


(10) 


until  a  convergence  criterion  is  reached.  The  following  two  conditions  are  required  of  the  SOS  distribution:  (i) 
ck\s  =  0  and  (ii)  ck+1  is  smooth.  To  fulfill  (ii),  a  regularization  term  similar  to  the  one  used  in  Eq.  (2)  was 
included.  The  step  size,  ek ,  can  be  determined  by  using  Armijo-Golstein  rule  or  by  simply  setting  ek  =  e.  The 
update  scheme  described  in  Eq.  (10)  takes  a  large  number  of  iterations  to  converge.  Therefore,  we  used  the 
limited-memory  Broydon-Fletcher  Goldfarb-Shanno  (L-BFGS)  method  to  solve  for  this  nonlinear  optimization 
problem. 


3.  RESULTS 

We  performed  image  reconstruction  simulations  for  two  different  numerical  phantoms  using  a  ray-tracing  method. 
In  one  of  these  breast  phantoms,  the  SOS  value  for  the  subcutaneous  fat  was  chosen  to  be  1375  m/s  and  the 
background  SOS  value  was  fixed  at  1480  m/s  to  model  a  large  variation  in  the  SOS  distribution.  In  the  second 
2D  breast  phantom,  SOS  value  for  the  subcutaneous  fat  was  chosen  to  be  1475  m/s  and  the  background  SOS 
value  was  chosen  to  bel500  m/s  to  yield  a  sample  with  small  contrast  in  the  SOS  distribution.  We  used  128 
sources  and  128  receivers  and  calculated  the  measured  TOF  data  using  the  bent-ray  method.  In  this  case,  if  the 
output  of  Eq.  (4)  is  close  to  the  exact  gradient  of  the  objective  function,  it  is  expected  to  recontrust  the  original 
SOS  distribution.  The  reconstructed  image  will  be  degraded  for  the  larger  variation  in  SOS  distribution  when 
the  linearized  approximation  is  used  to  calculate  V(c)  given  in  Eq.  (4). 

It  can  clearly  be  seen  from  the  Fig.  (1)  that  image  quality  is  degraded  when  the  contrast  in  the  SOS  distribu¬ 
tion  becomes  larger.  For  high  contrast  situations,  the  reconstructed  image  is  blurred  and  we  cannot  quantitatively 
recover  the  SOS  distribution.  In  principle,  both  images  will  look  like  the  numerical  phantom  if  the  gradient  in 
Eq.  (4)  is  evaluated  accurately.  We  will  present  results  from  the  adjoint-state  to  address  this  issue. 


4.  NEW  BREAKTHROUGH  WORK 

Our  formulation  for  the  ray-tracing  method  and  adjoint  state  method  will  be  applied  for  the  reconstruction  of 
speed  of  sound  distributions.  Our  focus  is  to  compare  these  methods  specifically  for  breast  imaging.  To  achieve 
this,  we  implemented  a  numerical  ring  scanner  geometry  similar  to  [N.  Duric  et.  at,  Med.  Phys.  34(2),  773 
(2007)].  We  will  present  the  comparison  for  a  numerical  breast  phantom  comprised  of  a  distribution  of  malignant, 
cystic,  and  fatty  masses  in  the  region  of  glandular  tissue  inside  a  ring  of  subcutaneous  fat  tissue.  Moreover,  we 
will  also  present  the  comparison  study  for  a  numerical  phantom  of  SOS  distribution  synthesized  from  a  slice 
through  an  MRI  of  a  breast.  The  comparison  between  these  method  will  be  quantified  by  assessing  the  quality 
of  images  using  some  quantitative  physical  tests  e.g.  least-square  error,  resolution. 


Figure  1:  Comparison  of  the  SOS  image  reconstruction  for  two  different  phantoms.  The  true  phantoms  are 
shown  in  the  left  column.  The  top  phantom  has  a  larger  variation  in  the  speed-of-sound  values.  The  middle 
column  of  images  shows  the  images  reconstructed  by  use  of  the  bent-ray  model,  with  the  corresponding  image 
profiles  displayed  in  the  right  column. 


5.  CONCLUSION 

This  study  demonstrates  the  use  of  the  adjoint-state  method  for  ultrasound  imaging.  The  ray-tracing  method 
has  primarily  been  used  in  ultrasound  computed  tomography,  but  is  cumbersome  and  prevents  the  accurate 
calculation  of  the  gradient  of  the  cost  function.  By  use  of  the  adjoint  method,  the  gradient  is  calculated  accurately 
and  improved  reconstructions  of  an  object’s  speed  of  sound  distribution  can  be  obtained. 
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ABSTRACT 

Optoacoustic  tomography  (OAT),  also  known  as  photoacoustic  computed 
tomography,  has  found  many  biomedical  applications.  Because  they  can 
model  complicated  imaging  physics,  compensate  for  imperfect  data 
acquisition  systems,  and  exploit  prior  information  regarding  the  object, 
iterative  image  reconstruction  algorithms,  in  general,  produce  higher  quality 
images  than  do  analytical  image  reconstruction  algorithms.  However,  three- 
dimensional  (3D)  iterative  image  reconstruction  is  computationally 
burdensome.  Even  with  graphics  processing  unit  (GPU)-accelerated 
implementations,  to  our  knowledge,  it  still  takes  at  least  five  hours  to 
reconstruct  the  3D  volume  of  a  whole-body  mouse.  This  computational 
burden  greatly  hinders  the  application  of  advanced  image  reconstruction 
algorithms  to  applications  with  a  large  field-of-view  (FOV),  such  as  breast 
imaging. 

In  this  study,  an  improved  GPU-based  implementation  of  a  numerical 
imaging  model  and  its  adjoint  have  been  developed  for  use  with  general 
gradient-based  iterative  image  reconstruction  algorithms.  Particularly,  two 
types  of  computation-reduced  discretization  methods  are  employed;  a 
parallel  fast  Fourier  transform  (FFT)  algorithm  is  employed  to  accelerate  the 
calculation  of  the  temporal  convolution  with  ultrasonic  transducer  responses; 
and  a  volume-reduction  method  is  proposed  to  reduce  the  computation  for 
applications  with  irregular  FOV.  Both  computer-simulation  and 
experimental  studies  are  conducted  to  investigate  the  efficiency  and 
accuracy  of  the  proposed  implementation.  The  results  suggest  that  the 
proposed  implementation  is  more  than  five  times  faster  than  previous 
implementations.  Using  the  proposed  implementation,  a  3D  whole-body 
mouse  image  can  be  reconstructed  in  less  than  one  hour.  The  developed 


algorithm  is  also  evaluated  for  3D  OAT  breast  imaging  with  sub  millimeter 
resolution. 

Keywords:  Optoacoustic  tomography,  iterative  image  reconstruction,  GPU 
acceleration,  unmatched  backprojection 
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ABSTRACT 

Ultrasound  computed  tomography  (USCT)  holds  great  promise  for  improving  the  detection  and  management  of 
breast  cancer.  Because  they  are  based  on  the  acoustic  wave  equation,  waveform  inversion-based  reconstruction 
methods  can  produce  images  that  possess  improved  spatial  resolution  properties  over  those  produced  by  ray- 
based  methods.  However,  waveform  inversion  methods  are  computationally  demanding  and  have  not  been  applied 
widely  in  USCT  breast  imaging.  A  computationally  efficient  numerical  wave  equation  solver  has  been  reported 
based  on  a  modified  Fresnel  propagation,  which  only  applies  to  USCT  systems  with  a  planar  incident  wave.  For 
breast  imaging  systems  with  a  spherical  incident  wave,  waveform  inversion-based  reconstruction  methods  remain 
computationally  challenging. 

In  this  work,  source  encoding  concepts  are  employed  to  develop  an  accelerated  USCT  reconstruction  method 
that  circumvents  the  large  computational  burden  of  conventional  waveform  inversion  methods.  This  method, 
referred  to  as  the  waveform  inversion  with  source  encoding  (WISE)  method,  encodes  the  measurement  data 
using  a  random  encoding  vector  and  determines  an  estimate  of  the  speed-of-sound  distribution  by  solving  a 
stochastic  optimization  problem  by  use  of  a  stochastic  gradient  descent  algorithm.  For  practical  applications,  a 
data-filling  strategy  is  proposed  to  mitigate  source  inferences  to  its  neighbor  receivers.  Computer-simulation  and 
experimental  phantom  studies  are  conducted  to  demonstrate  the  use  of  the  WISE  method.  Using  a  single  graphics 
processing  unit  card,  each  iteration  can  be  completed  within  25  seconds  for  a  128  x  128  mm2  reconstruction  region. 
The  results  suggest  that  the  WISE  method  maintains  the  high  spatial  resolution  of  waveform  inversion  methods 
while  significantly  reducing  the  computational  burden. 

1.  PURPOSE 

This  study  is  focused  on  the  image  reconstruction  of  breast  speed-of-sound  (SOS)  distribution  in  USCT.  The 
majority  of  USCT  image  reconstruction  methods  for  breast  imaging  investigated  to  date  have  been  based  on 
approximations  to  the  acoustic  wave  equation.1,2  A  relatively  popular  class  of  methods  is  based  on  geometrical 
acoustics.  They  are  commonly  referred  to  as  ‘ray-based’  methods.  Although  ray-based  methods  can  be  compu¬ 
tationally  efficient,  the  spatial  resoultion  of  the  images  they  produce  is  limited  due  to  the  fact  that  diffraction 
effects  are  not  modelled.3,4  This  is  undesirable  for  breast  imaging  applications,  in  which  the  ability  to  resolve 
fine  features,  e.g.,  tumor  spiculations,  is  important  for  distinguishing  healthy  from  diseased  tissues. 

USCT  reconstruction  methods  based  on  the  acoustic  wave  equation,  also  known  as  full-wave  inverse  scattering 
or  waveform  inversion  methods,  have  also  been  explored  for  a  variety  of  applications  including  medical  imaging. 4-7 
Because  they  account  for  higher-order  diffraction  effects,  waveform  inversion  methods  can  produce  images  that 
possess  higher  spatial  resolution  properties  than  those  produced  by  ray-based  methods.4,5  However,  conventional 
waveform  inversion  methods  are  iterative  in  nature  and  require  the  wave  equation  to  be  solved  numerically  a  large 
number  of  times  at  each  iteration.  Consequently,  such  methods  can  be  extremely  computationally  burdensome. 
For  special  geometries,7  efficient  numerical  wave  equation  solvers  have  been  reported.  However,  apart  from 
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special  cases,  the  large  computational  burden  of  waveform  inversion  methods  has  hindered  their  widespread 
application. 

The  purpose  of  this  study  is  to  develop  an  algorithmically  accelerated  waveform  inversion  method  for  breast 
SOS  reconstruction.  Aided  by  a  graphics  processing  unit  (GPU)-accelerated  implementation,  the  developed 
method  will  maintain  the  high  spatial  resolution  of  standard  waveform  inversion  methods  with  a  signifcant 
reduction  in  computational  time. 


2.  METHODS 

A  conventional  waveform  inversion  method  seeks  the  solution  of 

M-l 

c  =  arg nun  —  ^  ||g™  -  Hcsm||2  +  /311(c),  (1) 

m— 0 

where  c  is  the  sought-after  object  to  be  reconstructed,  i.e,  SOS  distribution,  gm  denotes  the  measured  data 
vector,  sm  denotes  the  (known)  source  vector,  Hc  denotes  a  numerical  wave  equation  solver  (NWES)  that 
maps  the  known  source  vector  to  the  measured  data  vector,  and  77(c)  and  (3  denote  the  penalty  term  and  the 
regularization  parameter  respectively.  The  superscript  in  Hc  indicates  the  dependence  of  Hc  on  c.  Note  that  one 
USCT  measurement  involves  firing  a  sequence  of  acoustic  pulses  in  turn  and  recording  the  data  corresponding  to 
every  pulse.  Each  pulse-firing  and  data  recording  process  will  be  indexed  by  m  for  m  =  0, 1,  ■  •  •  ,  M  —  1.  Solving 
Eqn.  (1),  in  general,  requires  the  calculation  of  ^  o  ^c||gm— Hcsm||2,  where  Vc  denotes  the  gradient  operator 
with  respect  to  c.  The  gradient  in  each  summand  is  commonly  computed  by  use  of  an  adjoint  state  method,5 
which  requires  two  runs  of  the  NWES.  Repeating  the  gradient  calculation  for  all  sources  results  in  2 M  runs 
of  the  NWES  at  each  iteration.  This  computational  burden  largely  hinders  the  application  of  the  conventional 
waveform  inversion  methods  in  practice. 

In  this  study,  a  waveform  inversion  with  source  encoding  (WISE)  method  was  developed.  The  WISE  method 
employs  the  objective  function 


c  =  argininEw{i||gw  -  Hcsw||2}  + /371(c),  (2) 

where  Ew  denotes  the  expectation  operator  with  respect  to  the  random  source  encoding  vector  w  £  Mm  ,  and 
gw  and  sw  denote  the  w-encoded  data  and  source  vectors,  defined  as 


M—l  M-l 

gw  =  ^  [w]mgm,  and  sw=5>]msm,  (3) 

m— 0  m— 0 

respectively.  Equation  (2)  was  solved  by  use  of  a  stochastic  gradient  descent  algorithm.8  Because  the  stochastic 
gradient  descent  algorithm  calculated  the  gradient  of  only  one  realization  of  the  random  variable  i  ||gw^MHcsw||2 
at  each  iteration,  the  required  number  of  NWES  runs  per  iteration  was  reduced  from  2 M  to  2.  Although  it,  in 
general,  requires  more  algorithm  iterations  to  average  out  the  randomness  in  the  realizations,  the  WISE  method, 
as  demonstrated  later,  can  greatly  reduce  the  overall  number  of  NWES  runs.  Both  computer-simulation  and 
experimental  phantom  studies  were  conducted  to  demonstrate  the  use  of  the  WISE  method  for  breast  SOS 
reconstruction. 


3.  RESULTS 

The  images  reconstructed  from  the  computer-simulated  noise-free  data  by  use  of  the  WISE  method  after  199 
iterations  and  sequential  waveform  inversion  method  after  43  iterations  are  shown  in  Fig.  l-(a)  and  (b).  As 
expected,4,9  both  images  are  more  accurate  and  possess  higher  spatial  resolution  than  the  one  reconstructed  by 
use  of  the  bent-ray  reconstruction  algorithm  displayed  in  Fig.  l-(c).  The  images  shown  in  Fig.  l-(a)  and  -(b) 
possess  similar  accuracies  as  measured  by  their  Euclidean  distances  from  the  SOS  phantom  vector  c,  namely 
0.07%  of  || c ||  for  the  former  and  0.08%  of  ||c||  for  the  latter.  However,  the  reconstruction  of  Fig.  l-(a)  required 


only  about  1.7%  of  the  computational  time  required  to  reconstruct  Fig.  l-(b),  namely,  1.4  hours  for  the  former 
and  81.4  hours  for  the  latter  respectively.  This  is  because  the  WISE  method  required  only  1018  NWES  runs, 
which  is  signficantly  less  than  the  58880  NWES  runs  required  by  the  sequential  waveform  inversion  method. 
With  a  similar  number  of  NWES  runs,  (e.g.,  1024),  one  can  only  complete  a  single  algorithm  iteration  by  use 
of  the  sequential  waveform  inversion  method.  The  corresponding  image,  shown  in  Fig.  l-(d),  lacks  quantitative 
accuracy  as  well  as  qualitative  value  for  identifying  features.  The  results  suggest  that  the  WISE  method  maintains 
the  advantages  of  the  sequential  waveform  inversion  method  while  significantly  reducing  the  computational  time. 


(a)  (b)  (c)  (d) 

Figure  1.  Images  reconstructed  by  use  of  (a)  the  WISE  method  after  the  199-th  iteration  (1,018  runs  of  NWES)  (b)  the 
sequential  waveform  inversion  algorithm  after  the  43-rd  iteration  (58,880  runs  of  the  NWES),  (c)  the  bent-ray  model- 
based  SOS  reconstruction  method,  and  (d)  the  sequential  waveform  inversion  algorithm  after  the  1-st  iteration  (1,024 
runs  of  the  NWES)  from  the  noise-free  non-attenuated  data.  The  grayscale  window  is  [1.46, 1.58]  mm//rs. 


The  images  reconstructed  from  the  experimentally-measured  data  are  shown  in  Fig.  2.  The  spatial  resolution 
of  the  image  reconstructed  by  use  of  the  WISE  method  is  significantly  higher  than  that  reconstructed  by  use 
of  the  bent-ray  model-based  method.  In  particular,  the  structures  labeled  ‘A’  and  ‘B’  possess  clearly-defined 
boundaries.  In  addition,  the  structure  labeled  ‘Cancer’  in  Fig.  2-(a)  is  almost  indistinguishable  in  the  image 
reconstructed  by  use  of  the  bent-ray  model-based  method  (see  Fig.  2-(b)).  The  improved  spatial  resolution  is 
expected  because  the  WISE  method  takes  into  account  the  high-order  diffraction  effects,  which  are  ignored  by 
the  bent-ray  method.4 


Figure  2.  (a)  Schematic  of  the  breast  phantom  employed  in  the  experimental  study.  Images  reconstructed  from  the 
experimentally  measured  phantom  data  by  use  of  (a)  the  bent-ray  model-based  SOS  reconstruction  method  and  (b)  the 
WISE  method  after  the  200-th  iteration.  The  grayscale  window  is  [1.49,  1.57]  mm//is. 


4.  NEW  OR  BREAKTHROUGH  WORK  TO  BE  PRESENTED 

Source  encoding  concepts  are  demonstrated  in  breast  USCT  experimental  studies  for  the  first  time.  Unlike 
previously  studied  waveform  inversion  methods  that  were  based  on  the  Helmholtz  equation,  the  WISE  method 


is  formulated  by  use  of  the  time-domain  acoustic  wave  equation.  A  GPU-accelerated  NWES  is  developed  that 
can  compute  1800  time  samples,  on  a  1024  x  1024  spatial  grid,  in  5  seconds.  In  addition,  a  data-filling  strategy 
is  proposed  to  mitigate  the  inference  of  the  source  with  its  neighboring  receivers  for  practical  applications. 

5.  CONCLUSION 

It  is  known  that  waveform  inversion-based  reconstruction  methods  can  produce  SOS  images  that  possess  im¬ 
proved  spatial  resolution  properties  over  those  produced  by  ray-based  methods.  However,  waveform  inversion 
methods  are  computationally  demanding  and  have  not  been  applied  widely  in  USCT  breast  imaging.  In  this 
work,  based  on  the  time-domain  wave  equation  and  motivated  by  recent  mathematical  results  in  the  geophysics 
literature,  the  WISE  method  was  developed  that  circumvents  the  large  computational  burden  of  conventional 
waveform  inversion  methods.  This  method  encodes  the  measurement  data  using  a  random  encoding  vector  and 
determines  an  estimate  of  the  speed-of-sound  distribution  by  solving  a  stochastic  optimization  problem  by  use 
of  a  stochastic  gradient  descent  algorithm.  With  our  current  GPU-based  implementation,  the  computation  time 
was  reduced  from  weeks  to  hours.  The  WISE  method  was  systematically  investigated  in  computer-simulation 
and  experimental  studies  involving  a  breast  phantom.  The  results  suggest  that  the  method  holds  value  for  USCT 
breast  imaging  applications  in  a  practical  setting. 

6.  DISCLOSURE 

This  work  is  original.  Parts  of  this  work  have  been  submitted  to  IEEE  Transactions  on  Ultrasounics,  Ferro- 
electrics  and  Frequency  Control  and  are  under  review. 
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Abstract 

Ultrasound  computed  tomography  (USCT)  holds  great  promise  for  improving  the  detection  and 
management  of  breast  cancer.  Because  they  are  based  on  the  acoustic  wave  equation,  waveform  inversion- 
based  reconstruction  methods  can  produce  images  that  possess  improved  spatial  resolution  properties 
over  those  produced  by  ray-based  methods.  However,  waveform  inversion  methods  are  computationally 
demanding  and  have  not  been  applied  widely  in  USCT  breast  imaging.  In  this  work,  source  encoding 
concepts  are  employed  to  develop  an  accelerated  USCT  reconstruction  method  that  circumvents  the 
large  computational  burden  of  conventional  waveform  inversion  methods.  This  method,  referred  to  as 
the  waveform  inversion  with  source  encoding  (WISE)  method,  encodes  the  measurement  data  using 
a  random  encoding  vector  and  determines  an  estimate  of  the  sound  speed  distribution  by  solving  a 
stochastic  optimization  problem  by  use  of  a  stochastic  gradient  descent  algorithm.  Both  computer- 
simulation  and  experimental  phantom  studies  are  conducted  to  demonstrate  the  use  of  the  WISE  method. 

The  results  suggest  that  the  WISE  method  maintains  the  high  spatial  resolution  of  waveform  inversion 
methods  while  significantly  reducing  the  computational  burden. 
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I.  Introduction 

After  decades  of  research  [1] — [4],  advancements  in  hardware  and  computing  technologies 
are  now  facilitating  the  clinical  translation  of  ultrasound  computed  tomography  (USCT)  for 
breast  imaging  applications  [5]— [9].  USCT  holds  great  potential  for  improving  the  detection  and 
management  of  breast  cancer  since  it  provides  novel  acoustic  tissue  contrasts,  is  radiation-  and 
breast-compression-free,  and  is  relatively  inexpensive.  [10],  [11].  Several  studies  have  reported 
the  feasibility  of  USCT  for  characterizing  breast  tissues  [4]— [7],  [11],  [12].  Although  some  USCT 
systems  are  capable  of  generating  three  images  that  depict  the  breast’s  acoustic  reflectivity, 
acoustic  attenuation,  and  sound  speed  distributions,  this  study  will  focus  on  the  reconstruction 
of  the  sound  speed  distribution. 

A  variety  of  USCT  imaging  systems  have  been  developed  for  breast  sound  speed  imaging 
[6],  [8],  [11],  [13]— [16].  In  a  typical  USCT  experiment,  acoustic  pulses  that  are  generated  by 
different  transducers  are  employed,  in  turn,  to  insonify  the  breast.  The  resulting  wavefieled  data 
are  measured  by  an  array  of  ultrasonic  transducers  that  are  located  outside  of  the  breast.  Here 
and  throughout  the  manuscript,  a  transducer  that  produces  an  acoustic  pulse  will  be  referred 
to  as  an  emitter;  the  transducers  that  receive  the  resulting  wavefield  data  will  be  referred  to  as 
receivers.  From  the  collection  of  recorded  wavefield  data,  an  image  reconstruction  method  is 
utilized  to  estimate  the  sound  speed  distribution  within  the  breast  [6],  [8],  [11]. 

The  majority  of  USCT  image  reconstruction  methods  for  breast  imaging  investigated  to  date 
have  been  based  on  approximations  to  the  acoustic  wave  equation  [13],  [17]-[24].  A  relatively 
popular  class  of  methods  is  based  on  geometrical  acoustics,  and  are  commonly  referred  to  as  ‘ray- 
based’  methods.  These  methods  involve  two  steps.  First,  time-of-flight  (TOF)  data  corresponding 
to  each  emitter-receiver  pair  are  estimated  [25].  Under  a  geometrical  acoustics  approximation, 
the  TOF  data  are  related  to  the  sound  speed  distribution  via  an  integral  geometry,  or  ray-based, 
imaging  model  [17],  [26].  Second,  by  use  of  the  measured  TOF  data  and  the  ray -based  imaging 
model,  a  reconstruction  algorithm  is  employed  to  estimate  the  sound  speed  distribution.  Although 
ray-based  methods  can  be  computationally  efficient,  the  spatial  resoultion  of  the  images  they 
produce  is  limited  due  to  the  fact  that  diffraction  effects  are  not  modelled  [23],  [27],  This  is 
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undesirable  for  breast  imaging  applications,  in  which  the  ability  to  resolve  fine  features,  e.g., 
tumor  spiculations,  is  important  for  distinguishing  healthy  from  diseases  tissues. 

USCT  reconstruction  methods  based  on  the  acoustic  wave  equation,  also  known  as  full-wave 
inverse  scattering  or  waveform  inversion  methods,  have  also  been  explored  for  a  variety  of 
applications  including  medical  imaging  [13],  [22],  [23],  [28]  and  geophysics  [29]— [3 1].  Because 
they  account  for  higher-order  diffraction  effects,  waveform  inversion  methods  can  produce  images 
that  possess  higher  spatial  resolution  than  those  produced  by  ray-based  methods  [23],  [28]. 
However,  conventional  waveform  inversion  methods  are  iterative  in  nature  and  require  the  wave 
equation  to  be  solved  numerically  a  large  number  of  times  at  each  iteration.  Consequently, 
such  methods  can  be  extremely  computationally  burdensome.  For  special  geometries  [13],  [32], 
efficient  numerical  wave  equation  solvers  have  been  reported.  However,  apart  from  special  cases, 
the  large  computational  burden  of  waveform  inversion  methods  has  hindered  their  widespread 
application. 

A  natural  way  to  reduce  the  computational  complexity  of  the  reconstruction  problem  is  to 
reformulate  it  in  a  way  that  permits  a  reduction  in  the  number  of  times  the  wave  equation 
needs  to  be  solved.  In  the  geophysics  literature,  source  encoding  methods  have  been  proposed 
to  achieve  this  [29]— [3 1  ] .  When  source  encoding  is  employed,  at  each  iteration  of  a  prescribed 
reconstruction  algorithm,  all  of  the  acoustic  sources  produced  by  the  emitters  are  combined 
(or  ‘encoded’)  by  use  of  a  random  encoding  vector;  So  are  the  measured  wavefield  data.  As  a 
result,  the  wave  equation  may  need  to  be  solved  as  few  as  twice  at  each  algorithm  iteration. 
In  conventional  waveform  inversion  methods,  this  number  would  be  equal  to  twice  the  number 
of  emitters  employed.  Although  conventional  waveform  inversion  methods  may  require  fewer 
algorithm  iterations  to  obtain  a  specified  image  accuracy  compared  to  source  encoded  methods, 
as  demonstrated  later,  the  latter  can  greatly  reduce  the  overall  number  of  times  the  wave  equation 
needs  to  be  solved. 

In  this  study,  a  waveform  inversion  with  source  encoding  (WISE)  method  for  USCT  sound 
speed  reconstruction  is  developed  and  investigated  for  breast  imaging  with  a  circular  transducer 
array.  The  WISE  method  determines  an  estimate  of  the  SOS  distribution  by  solving  a  stochastic 
optimization  problem  by  use  of  a  stochastic  gradient  descent  algorithm  [30],  [33].  Unlike  previ¬ 
ously  studied  waveform  inversion  methods  that  were  based  on  the  Helmholtz  equation  [22],  [23], 
the  WISE  method  is  formulated  by  use  of  the  time-domain  acoustic  wave  equation  [34]— [36] 
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and  utilizes  broad-band  measurements.  The  wave  equation  is  solved  by  use  of  a  computationally 
efficient  k-space  numerical  wave  equation  solver  that  is  accelerated  using  graphics  processing 
units  (GPUs).  In  order  to  mitigate  the  interference  of  the  emitter  on  its  neighboring  receivers,  a 
heuristic  data  replacement  strategy  is  proposed.  The  method  is  validated  in  computer-simulation 
studies  that  include  modelling  errors  and  other  physical  factors.  The  practical  applicability  of 
the  method  is  further  demonstrated  in  studies  involving  experimental  breast  phantom  data. 

The  remainder  of  the  paper  is  organized  as  follows.  In  Section  II,  USCT  imaging  models  in 
their  continuous  and  discrete  forms  are  reviewed.  A  conventional  waveform  inversion  method  and 
the  WISE  method  for  sound  speed  reconstruction  are  formulated  in  Section  III.  The  computer- 
simulation  studies  and  corresponding  numerical  results  are  presented  in  Sections  IV  and  V, 
respectively.  In  Section  VI,  the  WISE  method  is  further  validated  in  experimental  breast  phantom 
studies.  Finally,  the  paper  concludes  with  a  discussion  in  Section  VII. 

II.  Background:  USCT  imaging  models 

In  this  section,  imaging  models  that  provide  the  basis  for  image  reconstruction  in  waveform 
inversion-based  USCT  are  reviewed  in  their  continuous  and  discrete  forms. 


A.  USCT  imaging  model  in  its  continuous  form 


Although  a  digital  imaging  system  is  properly  described  as  a  continuous-to-discrete  (C-D) 
maping  (See  Chapter  7  in  [37]),  for  simplicity,  a  USCT  imaging  system  is  initially  described  in 
its  continuous  form  below. 

In  USCT  breast  imaging,  a  sequence  of  acoustic  pulses  is  transmitted  through  the  breast. 
We  denote  each  acoustic  pulse  by  sm(r,t)  £  Lr(M3  x  [0,  oo)),  where  each  pulse  is  indexed  by 
an  integer  m  for  m  =  0, 1,  •  •  •  ,  M  —  1  with  M  denoting  the  total  number  of  acoustic  pulses. 
Although  it  is  spatially  localized  at  the  emitter  location,  each  source  can  be  expressed  as  a 
function  of  space  and  time.  When  the  m-th  pulse  propagates  through  the  breast,  it  generates  a 
pressure  wavefield  distribution  denoted  by  pm{ r,  t )  e  L2(R3  x  [0,  oo)).  If  acoustic  absorption  and 
mass  density  variations  are  negligible,  pm(r,t)  in  an  unbounded  medium  satisfies  the  acoustic 
wave  equation  [381: 


W2pm{r,t) 


1  d 2 

c2(r)  dt 2 


Pm(T,t) 


-47rsm(r,t), 


(1) 
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where  c(r)  is  the  sought-after  sound  speed  distribution.  Equation  (1)  can  be  expressed  in  operator 
form  as 

Pm(r,t)  =ncsm(r,t),  (2) 

where  the  linear  operator  Hc  :  L2(M3  x  [0,  oo))  L2(M3  x  [0,  oo))  denotes  the  action  of  the  wave 
equation  and  is  independent  of  the  index  of  m.  The  superscript  ‘c’  indicates  the  dependence  of 
TLC  on  c(r). 

Consider  that  pm{ r,  t)  is  recorded  outside  of  the  object  for  r  E  Om  and  t  E  [0,  T],  where  Qm  C 
M3  denotes  a  continuous  measurement  aperture.  In  this  case,  when  discrete  sampling  effects  are 
not  considered,  the  USCT  imaging  model  can  be  described  as  a  continuous-to-continuous  (C-C) 
mapping  as: 

gm(r,t)  =  Mmncsm(r,t),  for  m  =  0, 1,  ■  ■  ■  ,  M  —  1,  (3) 

where  gm(r,t)  E  ILr  ( Qrn  x  [0,  T])  denotes  the  measured  data  function,  and  the  operator  A4rn  is 
the  restriction  of  'Hc  to  ilm  x  [0,  T\.  Introducing  the  m-dependent  operator  M.m  allows  Eqn.  (3) 
to  describe  USCT  imaging  systems  in  which  the  measurement  aperture  can  vary  with  emitter 
location.  Here  and  throughout  the  manuscript,  we  will  refer  to  the  process  of  firing  one  acoustic 
pulse  and  acquiring  the  corresponding  wavefield  data  as  one  data  acquisition  indexed  by  m.  The 
USCT  reconstruction  problem  in  its  continuous  form  is  to  estimate  the  sound  speed  distribution 
c(r)  by  use  of  Eqn.  (3)  and  the  data  functions  {gm(r, 

B.  USCT  imaging  model  in  its  discrete  forms 

A  digital  imaging  system  is  accurately  described  by  a  continuous-to-discrete  (C-D)  imaging 
model,  which  is  typically  approximated  in  practice  by  a  discrete-to-discrete  (D-D)  imaging  model 
to  facilitate  the  application  of  iterative  image  reconstruction  algorithms.  A  C-D  description  of  the 
USCT  imaging  system  is  provided  in  Appendix  A.  Below,  a  D-D  imaging  model  for  waveform 
inversion-based  sound  speed  reconstruction  in  USCT  are  presented.  The  D-D  imaging  model 
will  be  employed  subsequently  in  the  development  of  the  WISE  method  in  Section  III. 

Construction  of  a  D-D  USCT  imaging  model  requires  the  introduction  of  a  finite-dimensional 
approximate  representations  of  the  functions  c(r)  and  sm(r,t),  which  will  be  denoted  by  the 
vectors  c  E  M;V  and  sm  E  RNL.  Here,  N  and  L  denote  the  numbers  of  spatial  and  time  samples 
employed  for  wave  propagation  calculation  respectively.  In  waveform-based  USCT,  the  way  in 
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which  c(r)  and  sm(r,  t)  are  discretized  to  form  c  and  sm  is  dictated  by  the  numerical  method 
employed  to  solve  the  acoustic  wave  equation,  referred  to  as  a  numerical  solver.  In  this  study,  we 
employ  a  numerical  solver  based  on  a  pseudospectral  k-space  method  [34]— [36].  Accordingly, 
c(r)  and  sm(r,t)  are  sampled  on  Cartesian  grid  points  as 

[c]„  =  c(rn),  and  [sm\nL+i  =  sm(rn,  l A*),  for  (4) 

where  A*  denotes  the  temporal  sampling  interval  and  rn  denotes  the  location  of  the  n-th  point. 

For  a  given  c  and  sm,  the  pseudospectral  numerical  solver  can  be  described  in  operator  form 
as 

Pm  =  HCsm,  (5) 

where  the  matrix  Hc  is  of  dimension  NL  x  NL  and  represents  a  discrete  approximation  of  the 
wave  operator  LLC  defined  in  Eqn.  (2),  and  the  vector  p’’n  represents  the  estimated  pressure  data 
at  the  grid  point  locations  and  has  the  same  dimension  as  sm.  The  superscript  ‘a’  indicates  that 
these  values  are  approximate,  i.e.,  [p ^]nL+/  ~  /A*).  We  refer  the  readers  to  [34]— [36]  for 

additional  details  regarding  the  pseudospectral  numerical  solver. 

Because  the  pseudospectral  numerical  solver  yields  pressure  data  distributed  over  the  whole 
Cartesian  grid,  a  sampling  matrix  Mm  is  introduced  to  model  the  USCT  data  acquisition  process 
as 

grn  =  MmPm  =  MmHCSm,  (6) 


where  the  NTecL  x  NL  sampling  matrix  Mm  extracts  the  pressure  data  corresponding  to  the 
receiver  locations  on  the  measurement  aperture  Qm  with  Nrec  denoting  the  number  of  receivers, 
and  is  the  predicted  data  vector  that  approximates  the  true  measurements.  When  the  receiver 
and  grid  point  locations  do  not  coincide,  interpolation  methods  are  required.  As  an  example, 
when  a  nearest-neighbor  interpolation  method  is  employed,  the  elements  of  Mm  are  defined  as 


[Mm] 


niecL+l,nL+l 


1,  for  n  =  Zm(nrec), 
0,  otherwise, 


(7) 


where  \M.m]n^L+i,nL+i  denotes  the  element  of  Mm  at  the  ( nrecL  +  Z)-th  row  and  the  ( nL  +  Z)-th 
column,  and  lrn (nrec)  denotes  the  index  of  the  grid  point  that  is  closest  to  r (m,  nrec).  Here, 
r(m,  nrec)  denotes  the  location  of  the  nrec-th  receiver  at  the  m-th  data  acquisition.  In  summary, 
Eqn.  (6)  represents  the  D-D  imaging  model  that  will  be  employed  in  the  remainder  of  this  study. 
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Note  that  because  of  the  dependence  of  Mm  on  m,  a  varying  detection  geometry  among  data 
acquisitions  can  be  described  by  use  of  this  model. 


III.  Waveform  inversion  with  source  encoding  for  USCT 


A.  Sequential  waveform  inversion  in  its  discrete  form 


A  conventional  waveform  inversion  method  that  does  not  utilize  source  encoding  will  be 
employed  as  a  reference  for  the  developed  WISE  method  and  is  briefly  described  below.  Like 
other  conventional  approaches,  this  method  sequentially  processes  the  data  acquisitions  gm  for 
m  —  0, 1,  •  •  •  ,  M  —  1  at  each  iteration  of  the  associated  algorithm.  As  such,  we  will  refer  to  the 
conventional  method  as  a  sequential  waveform  inversion  method. 

A  sequential  waveform  inversion  method  can  be  formulated  as  a  non-linear  numerical  opti¬ 
mization  problem: 

c  =  argmin{Jr(c)  +  /57 £(c)},  (8) 

C 


where  ^(c),  7Z(c),  and  /3  denote  the  data  fidelity  term,  the  penalty  term,  and  the  regularization 
parameter,  respectively.  The  data  fidelty  term  J-(c)  is  defined  as  a  sum  of  squared  f2-norms  of 
the  data  residuals  corresponding  to  all  data  acquisitions  as: 

1  M-l 

F(c)  =  -  || Sm  -  MmHcsm||2,  (9) 

m= 0 

where  gm  6  L  denotes  the  measured  data  vector  at  the  m-th  data  acquisition.  The  choice 
of  the  penalty  term  will  be  addressed  in  Section  IV. 

The  gradient  of  T-'(c)  with  respect  to  c,  denoted  by  J,  will  be  computed  by  discretizing  an 
expression  for  the  Frechet  derivative  that  is  derived  assuming  a  continuous  form  of  Eqn.  (9). 
The  Frechet  derivative  is  described  in  Appendix  B.  Namely,  the  gradient  is  approximated  as 


M — 1  ivi — 1  J-j — z  Tal  nTal  iT 

r-ri  -  Vn  i  _  1  [Pm}nL+i-i  ~  2[pm]nL+i  +  [p 

[jjn  —  \n  ~  rc13  Z— /  /  AQmlnL+jL-l) 


M-l  L- 2 


>m. 


(10) 


m= 0 


m= 0  1=1 


where  Jm  denotes  the  gradient  of  |||gm  —  MmHcsm||2  with  respect  to  c  and  the  vector  q)), 
contains  samples  that  approximate  adjoint  wavefield  qm(r,  t )  that  satisfies  Eqn.  (33)  in  Appendix 
B.  By  use  of  the  pseudospectral  numerical  solver,  can  be  calculated  by 


qa 


mi 


(ii) 
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where 


r  !  _J  fe*&2WwL+(W).  if  n  G  Nm, 
[' T  m\nL+l  ^ 

0,  otherwise 


(12) 


Here,  Nm  =  {n  :  Xm(nrec),  nrec  =  0, 1,  •  •  •  ,  Nrec  —  1},  and  Im1  denotes  the  inverse  mapping  of 


r 

-L,m  • 


Given  the  explicit  form  of  J  in  Eqn.  (10),  a  variety  of  optimization  algorithms  can  be  employed 
to  solve  Eqn.  (8)  [39].  When  a  gradient  descent  algorithm  is  employed,  the  sequential  waveform 
inversion  method  is  given  by  Algorithm  1,  where  in  Line-10,  JK  denotes  the  gradient  of  72(c) 
with  respect  to  c. 


Algorithm  1  Gradient  descent-based  sequential  waveform  inversion. 
Input:  {g™},  { sm } ,  c® 

Output:  c 

l:  k  i —  0  \k  is  the  number  of  algorithm  iteration.} 

2:  while  stopping  criterion  is  not  satisfied  do 

3:  k  i —  k  -\-  1 

4:  J  <r~  0 

5:  for  m  :  =  0  to  M  —  1  do 

6:  Hcsm  {m  is  the  index  of  the  emitter.} 

7:  •<—  H cTm  (rm  is  calculated  via  Eqn.  (12).} 

8:  J  J  -f  Jm  { J,m  is  calculated  via  Eqn.  (10).} 

9:  end  for 

10:  J  <-  J  +  /3JR 

li:  Determine  step  size  A  via  a  line  search 

12:  C —  AJ 

13:  end  while 
14:  C  =  C ^ 


In  Algorithm  1,  Hc  is  the  most  computationally  burdensome  operator,  representing  one  run 
of  the  numerical  solver.  Note  that  it  appears  in  Lines-6,  -7,  and  -11.  Because  Lines-6  and  -7 
have  to  be  executed  M  times  to  process  all  of  the  data  acquisitions,  the  numerical  solver  has 
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to  be  executed  at  least  (2 M  +  1)  times  at  each  algorithm  iteration.  The  line  search  in  Line-11 
searches  for  a  step  size  along  the  direction  of  —  J  so  that  the  cost  function  is  reduced  by  use  of  a 
classic  trial-and-error  approach  [39].  Note  that,  in  general,  the  line  search  will  require  more  than 
one  application  of  Hc,  so  (2 M  +  1)  represents  a  lower  bound  on  the  total  number  of  numerical 
solver  runs  per  iteration. 

B.  Stochastic  optimization-based  waveform  inversion  with  source  encoding  (WISE) 

In  order  to  alleviate  the  large  computational  burden  presented  by  sequential  waveform  inversion 
methods  (e.g.,  Algorithm  1),  a  source  encoding  method  has  been  proposed  [22],  [29],  [401.  This 
method  has  been  formulated  as  a  stochastic  optimization  problem  and  solved  by  various  stochastic 
gradient-based  algorithms  [30],  [31].  In  this  section,  we  adapt  the  stochastic  optimization-based 
formulation  in  [30]  to  find  the  solution  of  Eqn.  (8). 

The  WISE  method  employs  the  same  cost  function  given  in  Eqn.  (8)  except  that  the  data 
fidelity  term  in  Eqn.  (9)  is  reformulated  as  the  expectation  of  a  random  quantity  as  [29]— [3 1], 


Algorithm  2  Waveform  inversion  with  source  encoding  (WISE)  algorithm. 

Input:  {g™},  {sm},  c(°) 

Output:  c 

l:  k  i —  0  \k  is  the  number  of  algorithm  iteration} 

2:  while  stopping  criterion  is  not  satisfied  do 

3:  k  <—  k  +  1 

4:  Draw  elements  of  w  from  independent  and  identical  Rademacher  distribution. 

5:  pw  •<—  Hcsw  {sw  is  calculated  via  Eqn.  (14).  } 

6:  qw  <—  H crw  {See  text  for  the  calculation  of  rw} 

7:  J  <(—  Jw  +  /5JR  {Jw  is  calculated  via  Eqn.  (16)} 

8:  Determine  step  size  A  by  use  of  line  search 

9:  <-  -  AJW 

10:  end  while 

11:  C  = 
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[33],  [40],  [41] 


J-s(c)  =  Ew{-||gw-MHcsw||2}, 


(13) 


where  Ew  denotes  the  expectation  operator  with  respect  to  the  random  source  encoding  vector 
w  e  RM,  M  =  Mm  is  the  sampling  matrix  that  is  assumed  to  be  identical  for  m  —  0, 1,  ■  ■  ■  ,  M— 
1,  and  gw  and  sw  denote  the  w-encoded  data  and  source  vectors,  defined  as 

M—l  M— 1 

gw  =  ^  [w]mgm,  and  sw  =  ^[w]msm,  (14) 

m=0  m=0 

respectively.  It  has  been  demonstrated  that  Eqns.  (9)  and  (13)  are  mathematically  equivalent 
when  w  possesses  a  zero  mean  and  an  identity  covariance  matrix  [30],  [33],  [41].  In  this  case, 
the  optimization  problem  whose  solution  specifies  the  sound  speed  estimate  can  be  re-expressed 
in  a  stochastic  framework  as 


c  =  argminEw{-||gw  —  MHcsw||2}  +  /371(c) 


(15) 


which  we  refer  to  as  the  waveform  inversion  with  source  encoding  (WISE)  method.  An  im¬ 
plementation  of  the  WISE  method  that  utilizes  the  stochastic  gradient  descent  algorithm  is 
summarized  in  Algorithm  2. 

In  Algorithm  2,  the  numerical  solver  needs  to  be  run  one  time  in  each  of  Lines  5  and  6.  In 
the  line  search  to  determine  the  step  size  in  Line  8,  the  numerical  solver  needs  to  be  run  at  least 
one  time,  but  in  general  will  require  a  small  number  of  additional  runs,  just  as  in  Algorithm  1. 
Accordingly,  the  lower  bound  on  the  number  of  required  numerical  solver  runs  per  iteration  is  3, 
as  opposed  to  (2M  +  1)  for  the  conventional  sequential  waveform  inversion  method  described  by 
Algorithm  1.  As  demonstrated  in  geophysics  applications  [29],  [31],  [40]  and  the  breast  imaging 
studies  below,  the  WISE  method  provides  a  substantial  reduction  in  reconstruction  times  over 
use  of  the  standard  sequential  waveform  inversion  method.  In  Line-7,  Jw  can  be  calculated 
analogously  to  Eqn.  (10)  as 


L—2 

Eh" 


\nL+(L—l) ' 


[pW]nL+Z-l 


2  [p  W]nL+l  +  [pW]ni+Z+l 


1=1 


At 


(16) 


where  pw  =  Hcsw  and  qw  =  Hctw.  Various  probablity  density  functions  have  been  proposed 
to  describe  the  random  vector  w  [29],  [31],  [40].  In  this  study,  we  employed  a  Rademacher 
distribution  as  suggested  by  [29],  in  which  case  each  element  of  w  had  a  50%  chance  of  being 
either  +1  or  —1. 
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IV.  Description  of  computer-simulation  studies 

Two-dimensional  computer-simulation  studies  were  conducted  to  validate  the  WISE  method 
for  breast  sound  speed  imaging  and  demonstrate  its  computational  advantage  over  the  standard 
sequential  waveform  inversion  method. 

A.  Measurement  geometry 

A  circular  measurement  geometry  was  chosen  to  emulate  a  previously  reported  USCT  breast 
imaging  system  [11],  [23],  [42].  As  sketched  in  Fig.  1,  256  ultrasonic  transducers  were  uniformly 
distributed  on  a  ring  of  radius  Rs  =  110  mm.  The  generation  of  one  USCT  data  set  consisted 
of  M  =  256  sequential  data  acquisitions.  In  each  data  acquisition,  one  emitter  produced  an 
acoustic  pulse.  The  acoustic  pulse  was  numerically  propagated  through  the  breast  phantom  and 
the  resulting  wavefield  data  were  recorded  by  all  transducers  in  the  array  as  described  below. 
Note  that  the  location  of  the  emitter  in  every  data  acquisition  was  different  from  those  in  other 
acquisitions,  while  the  locations  of  receivers  were  identical  for  all  acquisitions. 

B.  Numerical  breast  phantom 

A  numerical  breast  phantom  of  diameter  98  mm  was  employed.  The  phantom  was  composed 
of  8  structures  representing  adipose  tissues,  parenchymal  breast  tissues,  cysts,  benign  tumors,  and 
malignant  tumors,  as  shown  in  Fig.  2.  For  simplicity,  the  acoustic  attenuation  of  all  tissues  was 
described  by  a  power  law  with  a  fixed  exponent  y  —  1.5  [43].  The  corresponding  sound  speed 
values  and  the  attenuation  coefficients  are  listed  in  TABFE  I  [43]-[45].  Both  the  sound  speed 
and  the  absorption  coefficient  distributions  in  Fig.  2  were  sampled  on  a  uniform  Cartesian  grid 
with  spacing  As  =  0.25  mm.  The  finest  structure  (indexed  by  7  in  Fig.  2-(a))  was  of  diameter 
3.75  mm. 
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C.  Simulation  of  the  measurement  data 

1)  First-order  numerical  wave  equation  solver:  Acoustic  wave  propagation  in  acoustically 
absorbing  media  was  modeled  by  three  coupled  first-order  partial  differential  equations  [46]: 

7^u(M)  =  — ' Vp(r,f)  (17a) 

d 

— p(r,t)  =  — V  •  u(r,  t)  +  47t  /  dt's(r,t  )  (17b) 

ot  Jo 

p(r,t)  =  c2(r)[l  +  r(r)^(-V2)s//2_1  +r7(r)(-V2)(y+1}/2_1]p(r,f),  (17c) 

where  u(r,i),  p(r,t),  and  p(r)  denote  the  acoustic  particle  velocity,  the  acoustic  pressure,  and 
the  acoustic  density,  respectively.  The  functions  r(r)  and  77 (r)  describe  acoustic  absorption  and 
dispersion  during  the  wave  propagation  respectively  [46]: 

r(r)  =  -2«o(r)c0(r)y_1,  7?(r)  =  2a0(r)c0(r)y  tan(?n//2),  (18) 

where  o-ofr)  and  y  are  the  absorption  coefficient  and  the  power  law  exponent  respectively.  When 
the  medium  is  assumed  to  be  lossless,  i.e.,  a0(r)  =  0,  it  can  be  shown  that  Eqn.  (17)  is  equivalent 
to  Eqn.  (1). 

Based  on  Eqn.  (17),  a  pseudospectral  k-space  method  was  employed  to  simulate  acoustic 
pressure  data  [36],  [46].  The  numerical  scheme  will  be  referred  to  as  a  first-order  numerical 
solver.  The  first-order  numerical  solver  was  implemented  using  graphic  processing  units  (GPUs). 
The  calculation  domain  was  of  size  512  x  512  mm2,  sampled  on  a  2048  x  2048  uniform  Cartesian 
grid  of  spacing  As  =  0.25  mm.  A  nearest-neighbor  interpolation  was  employed  to  place  all 
transducers  on  the  grid  points.  On  a  platform  consisting  of  dual  quad-core  CPUs  with  a  3.30  GHz 
clock  speed,  64  gigabytes  (GB)  of  random- accessing  memory  (RAM),  and  a  single  NVIDIA  Tesla 
K20  GPU,  the  first-order  numerical  solver  took  108  seconds  to  complete  one  foward  simulation. 

2)  Acoustic  excitation  pulse:  The  excitation  pulse  employed  in  this  study  was  assumed  to  be 
spatially  localized  at  the  emitter  location  while  temporally  it  was  a  fc  =  0.8  MHz  sinusoidal 
function  tapered  by  a  Gaussian  kernel  with  standard  deviation  a  =  0.5  ps,  i.e., 

{exp  (  —  ~—y~T~  j  sin(27r  fct),  at  the  m-th  emitter  location 

(19) 

0,  otherwise. 

The  temporal  profile  and  the  amplitude  frequency  spectrum  of  the  excitation  pulse  are  plotted 
in  Fig.  3-(a)  and  -(b),  respectively.  The  excitation  pulse  contained  approximate  3  cycles. 
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3)  Generation  of  non-attenuated  and  attenuated  noise-free  data:  For  every  data  acquisition, 
indexed  by  m,  the  first-order  numerical  wave  equation  solver  was  run  for  3600  time  steps  with 
a  time  interval  A*  =  0.05/is  (corresponding  to  a  20  MHz  sampling  rate).  Downsampling  the 
recorded  data  by  taking  every  other  time  sample  resulted  in  a  data  vector,  i.e.,  gm  in  Eqn.  (9),  that 
was  effectively  sampled  at  10  MHz  and  was  of  dimensions  ML  with  M  =  256  and  L  =  1800. 
The  data  vector  at  the  0-th  data  acquisition,  g0,  is  displayed  as  a  2D  image  in  Fig.  4-(a).  This 
undersampling  procedure  was  introduced  to  avoid  inverse  crime  [47]  so  that  the  data  generation 
and  the  image  reconstruction  employed  different  numerical  discretization  schemes.  Repeating  the 
calculation  for  m  =  0, 1,  •  •  •  ,  255,  we  obtained  a  collection  {gm}  of  data  vectors  that  together 
represented  one  complete  data  set.  Utilizing  the  absorption  phantom  described  in  Section  IV-B, 
a  complete  attenuated  data  set  was  computed.  An  idealized,  non-attenuated,  data  set  was  also 
computed  by  setting  a0(r)  =  0. 

4)  Generation  of  incomplete  data:  An  incomplete  data  set  in  this  study  corresponds  to  one 
in  which  only  NTec  receivers  located  on  the  opposite  side  of  the  emitter  record  the  pressure 
wavefield,  with  Nrec  <  M.  Taking  the  0-th  data  acquisition  as  an  example  (see  Fig.  1),  only 
Nrec  =  100  receivers,  indexed  from  78  to  177,  record  the  wavefield,  while  other  receivers  record 
either  unreliable  or  no  measurements.  Incomplete  data  sets  formed  in  this  way  can  emulate  two 
practical  scenarios:  (1)  Signals  recored  by  receivers  near  the  emitter  are  unreliable  and  therefore 
discarded  [23];  and  (2)  An  arc-shaped  transducer  array  is  employed  that  rotates  with  the  emitter 
[14],  [15],  [48]. 

Specifically,  incomplete  data  sets  were  generated  as 


771=0, !,•••  ,M— 1 


for 


(20) 


nrec=0,l,---,iVrec  — O 


Jm(nrec)L+P 


where  g))'cpl  is  the  incomplete  m-th  data  acquisition,  which  is  of  dimensions  NTecL,  with  Nrec  < 
M.  The  index  map  Jm  :  {0, 1,  ■  •  •  ,  NTec  —  1}  i — M^od  is  defined  as 


mod  M, 


(21) 


where  ( m '  mod  M)  calculates  the  remainder  of  m!  divided  by  M,  and  the  index  set  MAKKi 
collects  indices  of  transducers  that  reliably  record  data  at  the  m-th  data  acquisition  and  is  defined 


as 


Msf,od  =  {k  mod  M\k  e  [m  +  (M  -  lVrec)/2,  m  +  (M  +  Nrec)/2)  j .  (22) 
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Here,  for  simplicity,  we  assume  both  M  and  Nrec  to  be  even  numbers.  In  this  study,  we 
empirically  set  Nrec  =  100  so  that  the  object  can  be  fully  covered  by  the  fan  region  as  shown 
in  Fig.  1. 

5)  Generation  of  noisy  data:  An  additive  Gaussian  white  noise  model  was  employed  to 
simulate  electronic  measurement  noise  as 

gm  =  gm  +  n,  (23) 

where  gm  and  n  are  the  noisy  data  vector  and  the  Gaussian  white  noise  vector,  respectively.  In 
this  study,  the  maximum  value  of  the  pressure  received  by  the  128-th  transducer  at  the  0-th  data 
acquisition  with  a  homogeneous  medium  (water  tank)  was  chosen  as  a  reference  signal  amplitude. 
The  noise  standard  deviation  was  set  to  be  5%  of  this  value.  An  example  of  a  simulated  noiseless 
and  noisy  data  acquisition  is  shown  Fig.  4. 

D.  Image  reconstruction 

1 )  Second-order  numerical  wave  equation  solver:  In  the  reconstruction  methods  described 
below,  the  action  of  the  operator  Hc  (Eqn.  (5))  was  computed  by  solving  Eqn.  (1)  by  use  of 
the  pseudospectral  k-space  numerical  solver  method.  This  was  implemented  using  GPUs.  The 
calculation  domain  was  of  size  512  x  512  mm2,  sampled  on  a  1024  x  1024  uniform  Cartesian 
grid  of  spacing  As  =  0.5  mm  for  reconstruction.  On  a  platform  consisting  of  dual  octa-core 
CPUs  with  a  2.00  GHz  clock  speed,  125  GB  RAM,  and  a  single  NVIDIA  Tesla  K20C  GPU, 
the  second-order  numerical  solver,  took  7  seconds  to  complete  one  forward  simulation. 

2)  Sequential  waveform  inversion:  To  serve  as  a  reference  for  the  WISE  method,  we  imple¬ 
mented  the  sequential  waveform  inversion  method  described  in  Algorithm  1.  No  penalty  term 
was  included  (/ 3  =  0)  because,  due  to  its  extreme  computational  burden,  we  only  investigated 
this  method  in  preliminary  studies  involving  noise-free  non-attenuated  data.  A  uniform  sound 
speed  distribution  was  employed  as  the  initial  guess.  We  assumed  that  the  background  sound 
speed  was  known  and  the  object  was  contained  in  a  square  region-of-interest  (ROI)  of  dimension 
128  x  128  mm2  (See  Fig.  1),  which  corresponded  to  256  x  256  pixels. 

3)  WISE  method:  We  implemented  the  WISE  method  according  to  Algorithm  2.  Two  types 
of  smoothness  penalties  were  employed  in  this  study:  a  quadratic  penalty  expressed  as 

77.Q(c)  =  'y  ^  'y  ^([cjjjv^-i-t  —  [c]jNx+i-i )  +  ([c }jNx+i  —  [c](j— l)JVj+i)  :  (24) 

j  i 
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where  Nx  and  Ny  denote  the  number  of  grid  points  along  the  ‘x’  and  ‘y’  directions  respectively, 
and  a  total  variation  (TV)  penalty,  defined  as  [49],  [50] 


^TV(c) =  5Z  5Z  \!e  +  (tcW+* -  Mitvx+i- 1)2  +  (Mjtvx+i  -  Ho-o^+i)2,  (25) 

3  i 

where  e  is  a  small  number  introduced  to  avoid  dividing  by  0  in  the  gradient  calculation.  In  this 
study,  we  empirically  selected  e  =  10-8.  The  value  was  fixed  because  we  observed  that  the  value 
of  e  had  a  minor  impact  on  the  reconstructed  images  compared  to  the  impact  of  (3.  In  addition, 
the  use  of  this  parameter  can  be  avoided  when  advanced  optimization  algorithms  are  employed 
[51].  As  in  the  sequential  waveform  inversion  case,  it  was  assumed  that  the  background  sound 
speed  was  known  and  the  object  was  contained  in  a  square  ROI  of  dimension  128  x  128  mm2  (See 
Fig.  1),  which  corresponded  to  256  x  256  pixels.  The  regularization  parameters  corresponding 
to  the  quadratic  penalty  and  the  TV  penalty  will  be  denoted  by  (3®  and  /!TV,  respectively. 
Optimal  regularization  parameter  values  depend  on  the  specific  medically  relavent  task  that  the 
reconstructed  images  are  used  for  [37].  Estimation  of  the  optimal  values  requires  a  systematic 
investigation  of  the  image  statistics,  which  is  out  of  the  scope  of  this  study.  Here,  we  only 
investigated  the  impacts  of  /)(-  and  /3TY  on  the  reconstructed  images  by  sweeping  their  values 
over  a  wide  range. 

4)  Reconstruction  from  incomplete  data:  Because  the  WISE  method  requires  Mm  to  be 
identical  for  all  m’s,  image  reconstruction  from  incomplete  data  remains  challenging  [30],  [33], 
[41].  In  this  study,  two  data  completion  strategies  were  investigated  [30],  [33],  [41]  to  synthesize 
a  complete  data  set,  from  which  the  WISE  method  could  be  effectively  applied. 

One  strategy  was  to  fill  the  missing  data  with  pressure  corresponding  to  a  homogeneous 
medium  as 


^combHi 


J  rarecL+Z 


[gmCPl] 


[g  m}m™L+h 

for  mrec  =  0, 1,  •  •  •  ,  M  —  1,  where  g,d  G  gm 


h  ~\Q)  ML  ^incpl  ^ 


if  mrec  G  M^od 

otherwise, 

^reci,  and  g“mbH  € 


(26) 


t>ML 


,  denote  the 


computer-simulated  (with  a  homogeneous  medium),  the  measured  incomplete,  and  the  combined 
complete  data  vectors  at  the  m-th  data  acquisition,  respectively.  The  mapping  Jmx  :  M^od  t-G 
{0, 1,  •  •  •  ,  Nrec  —  1}  denotes  the  inverse  operator  of  Jm  as 


27TKec)  = 


m 
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This  data  completion  strategy  is  based  on  the  assumption  that  the  back-scatter  from  breast  tissue 
in  an  appropriately  sound  speed-matched  water  bath  is  weak.  This  assumption  suggests  that  the 
missing  measurements  can  be  replaced  by  the  corresponding  pressure  data  that  would  have  been 
produced  in  the  absence  of  the  object. 

The  second,  more  crude,  data  completion  strategy  was  to  simply  fill  the  missing  data  with 
zeros,  i.e.. 


(28) 


where  g"mib0  denotes  the  data  completed  with  the  second  strategy. 

5)  Bent-ray  image  reconstruction:  A  bent-ray  method  was  also  employed  to  reconstruct  im¬ 
ages.  Details  regarding  the  time-of-flight  estimation  and  algorithm  implemenations  are  provided 
in  Appendix  C. 


V.  Computer-simulation  results 


A.  Images  reconstructed  from  idealized  data 

The  images  reconstructed  from  the  noise-free,  non-attenuated,  data  by  use  of  the  WISE  method 
with  199  iterations  and  the  sequential  waveform  inversion  method  with  43  iterations  are  shown 
in  Fig.  5-(a)  and  (b).  As  expected  [23],  [52],  both  images  are  more  accurate  and  possess  higher 
spatial  resolution  than  the  one  reconstructed  by  use  of  the  bent-ray  reconstruction  algorithm 
displayed  in  Fig.  5-(c).  Profiles  through  the  reconstructed  images  are  displayed  in  Fig.  6.  The 
images  shown  in  Fig.  5-(a)  and  -(b)  possess  similar  accuracies  as  measured  by  their  Euclidean 
distances  from  the  sound  speed  phantom  vector  c;  namely  0.07%  of  ||c||  for  the  former  and 
0.08%  of  || c||  for  the  latter.  However,  the  reconstruction  of  Fig.  5-(a)  required  only  about  1.7% 
of  the  computational  time  required  to  reconstruct  Fig.  5-(b);  namely,  1.4  hours  for  the  former 
and  81.4  hours  for  the  latter  respectively.  This  is  because  the  WISE  method  required  only 
1018  numerical  solver  runs  which  is  signficantly  less  than  the  57088  numerical  solver  runs 
required  by  the  sequential  waveform  inversion  method.  With  a  similar  number  of  numerical 
solver  runs,  (e.g.,  1024),  one  can  only  complete  a  single  algorithm  iteration  by  use  of  the 
sequential  waveform  inversion  method.  The  corresponding  image,  shown  in  Fig.  5-(d),  lacks 
quantitative  accuracy  as  well  as  qualitative  value  for  identifying  features.  The  results  suggest 
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that  the  WISE  method  maintains  the  advantages  of  the  sequential  waveform  inversion  method 
while  significantly  reducing  the  computational  time. 

B.  Convergence  of  the  WISE  method 

Images  reconstructed  from  noise-free,  non-attenuated,  data  by  use  of  the  WISE  method  contain 
radial  streak  artifacts  when  the  algorithm  iteration  number  is  less  than  100,  as  shown  in  Figs.  7- 
(a-c).  Profiles  through  these  images  are  displayed  in  8.  The  streaks  artifacts  are  likely  caused  by 
crosstalk  introduced  during  the  source  encoding  procedure  [31],  [40].  However,  these  artifacts 
are  effectively  mitigated  after  more  iterations  as  demonstrated  by  the  image  reconstructed  after 
the  199-th  iteration  in  Fig.  5-(a)  and  its  profile  in  Fig.  6.  The  quantitative  accuracy  of  the 
reconstructed  images  is  improved  with  more  iterations  as  shown  in  Fig.  8. 

Figure  9-(a)  reveals  that  the  WISE  method  requires  a  larger  number  of  algorithm  iterations 
than  does  the  sequential  waveform  inversion  method  to  achieve  the  same  reconstruction  accuracy, 
quantified  by  the  Euclidean  distance  of  the  reconstructed  image  in  percentage  of  ||c||.  Also,  the 
Euclidean  distance  of  the  reconstructed  images  by  use  of  the  WISE  method  appears  to  oscillate 
around  0.08%  of  ||c||  after  the  first  100  iterations  while  the  sequential  waveform  inversion  method 
can  achieve  a  higher  accuracy.  However,  as  shown  in  Fig.  5-(a),  the  image  reconstructed  by 
use  of  the  WISE  method  is  highly  accurate  after  the  199-th  iteration.  Moreover,  to  achieve 
the  same  accuracy,  the  amount  of  computation  for  the  WISE  method  is  about  two-order  of 
magnitude  smaller  than  that  for  the  sequential  waveform  inversion  method  as  suggested  by 
Fig.  9-(b).  This  is  because  of  the  significant  computation  reduction  per  iteration  when  the  WISE 
method  is  employed.  We  also  plotted  the  cost  function  value  against  the  number  of  iterations 
in  Fig.  9-(c).  Note  that  for  the  WISE  method,  the  cost  function  value  was  approximated  by  the 
current  realization  of  ^||gw  —  MHcsw||2.  These  plots  suggest  that,  in  this  particular  case,  the 
WISE  method  appears  to  approximately  converge  after  200  iterations.  For  example,  the  images 
reconstructed  after  199  (Fig.  5-(a))  and  250  (Fig.  7-(d))  iterations  are  nearly  identical. 

C.  Images  reconstructed  from  non-attenuated  data  containing  noise 

Images  reconstructed  by  use  of  the  WISE  method  with  a  quadratic  penalty  and  the  WISE 
method  with  a  TV  penalty  from  noisy,  non-attenutated,  data  are  presented  in  Fig.  10.  All 
images  were  obtained  after  1024  algorithm  iterations.  The  WISE  method  with  a  quadratic  penalty 
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effectively  mitigates  image  noise  as  shown  in  Figs.  lO-(a-c),  at  the  expense  of  image  resolution, 
as  expected.  Figure  10-(d)  shows  an  image  reconstructed  by  use  of  the  WISE  method  with  a 
TV  penalty.  The  image  appears  to  possess  a  similar  resolution  but  a  lower  noise  level  than 
the  image  in  Fig.  10-(b)  that  was  reconstructed  by  use  of  the  WISE  method  with  a  quadratic 
penalty. 

D.  Images  reconstructed  from  acoustically  attenuated  data 

Our  current  implementation  of  the  WISE  method  assumes  an  absorption-free  acoustic  medium. 
This  assumption  can  be  strongly  violated  in  practice.  In  order  to  investigate  the  robustness  of 
the  the  WISE  method  to  model  errors  associated  with  ignoring  medium  acoustic  absorption, 
we  applied  the  algorithm  to  the  acoustically  attenuated  data  that  were  produced  as  described 
in  Section  IV-C.  As  shown  in  Fig.  11,  when  the  medium  acoustic  absorption  is  considered,  the 
amplitude  of  the  measured  pressure  is  attenuated  by  approximately  a  factor  of  2.  The  wavefront 
(See  Fig.  ll-(a))  remains  very  similar  to  that  when  medium  absorption  is  ignored  (See  Fig.  4- 
(a)).  Medium  absorption  has  the  largest  impact  on  the  pressure  data  received  by  transducers 
located  opposite  the  emitter  as  shown  in  Fig.  ll-(b).  The  shape  of  the  pulse  profile  remains 
very  similar  as  shown  in  Fig.  ll-(c)  and  -(d),  suggesting  that  waveform  dispersion  may  be  less 
critical  than  amplitude  attenuation  in  image  reconstruction  for  this  phantom. 

Images  reconstructed  by  use  of  the  WISE  method  with  a  TV  penalty  from  noise-free  and 
noisy  attenuated  data  are  shown  in  Figs.  12-(a)  and  (b).  Image  profiles  are  shown  in  Fig.  12-(c). 
Although  these  images  contain  certain  artifacts  that  were  not  produced  in  the  idealized  data 
studies,  most  object  structures  remain  readily  identified.  These  results  suggest  that  the  WISE 
method  with  a  TV  penalty  can  tolerate  data  inconsistencies  associated  with  neglecting  acoustic 
attenuation  in  the  imaging  model,  at  least  to  a  certain  level  with  regards  to  feature  detection 
tasks. 

E.  Images  reconstructed  from  idealized  incomplete  data 

The  wavefront  of  the  noise-  and  attenuation-free  pressure  wavefield  when  the  object  is  absent 
(Fig.  13-(a))  appears  to  be  very  similar  to  that  when  the  object  is  present  (Fig.  4-(a)).  As 
expected,  the  largest  differences  are  seen  in  the  signals  received  by  the  transducers  located 
opposite  of  the  emitter,  as  shown  in  Fig.  13-(b).  As  seen  in  Fig.  13-(c),  the  time  traces  received 
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by  the  40-th  transducer  are  nearly  identical  when  object  is  present  and  absent.  This  is  because 
the  back-scattered  wavefield  is  weak  for  breast  imaging  applications.  These  results  suggest  the 
potential  efficacy  of  the  data  completion  strategy  of  filling  the  missing  data  with  the  pressure 
data  corresponding  to  a  water  bath. 

The  image  reconstructed  from  the  measurements  completed  with  pressure  data  corresponding 
to  a  water  bath  is  shown  in  Fig.l4-(a).  As  revealed  by  the  profile  in  Fig.l4-(c),  this  image  is 
highly  accurate.  Alternatively,  the  image  reconstructed  from  the  the  data  completed  with  zeros 
contains  strong  artifacts  as  shown  in  Fig.  14-(b).  These  results  suggest  that  the  WISE  method  can 
be  adapted  to  reconstruct  images  from  incomplete  data,  which  is  particularly  useful  for  emerging 
laser-induced  USCT  imaging  systems  [14]— [16]. 

VI.  Experimental  validation 

A.  Data  acquisition 

The  SoftVue  scanner  was  employed  in  the  experimental  study  [53].  The  scanner  contained 
a  ring-shaped  transducer  array  of  radius  110  mm.  2048  detecting  elements  were  uniformly 
distributed  on  the  ring.  Each  element  had  a  center  frequency  of  2.75  MHz  and  a  pitch  of 
0.34  mm.  Each  element  was  elevationally  focused  to  isolate  a  slice  of  3  mm  in  thickness.  The 
transducer  array  was  mounted  in  a  water  tank  and  could  be  translated  with  a  motorized  gantry 
in  the  vertical  direction.  We  refer  the  readers  to  [53]  for  more  details  about  the  USCT  imaging 
system. 

The  breast  phantom  was  built  by  Dr.  Ernie  Madsen  from  the  University  of  Wiscosin  and 
provides  tissue-equivalent  scanning  characteristics  of  highly  scattering,  predominantly  parenchy¬ 
mal  breast  tissue.  The  phantom  mimics  the  presence  of  benign  and  cancerous  masses  embedded 
in  glandular  tissue,  including  a  subcutaneous  fat  layer.  Figure  15  is  a  schematic  plot  of  one 
phantom  slice.  The  diameter  of  the  inclusions  is  approximately  12  mm.  Table  II  presents  the 
known  acoustic  properties  of  the  phantom. 

During  data  acquisition,  the  breast  phantom  was  placed  near  the  center  of  the  ring-shaped 
transducer  array  so  that  the  distance  between  the  phantom  and  each  transducer  was  approximately 
the  same.  While  scanning  each  slice,  every  other  transducer  element  sequentially  emits  fan  beam 
ultrasound  signals  towards  the  opposite  side  of  the  ring.  The  forward  scattered  and  backscattered 
ultrasound  signals  are  subsequently  recorded  by  the  same  transducer  elements.  The  received 
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waveform  was  sampled  at  a  rate  of  12  MHz.  The  1024  data  acquisitions  took  about  20  seconds 
in  total.  A  calibration  data  set  was  also  acquired  with  water  bath  only. 

B.  Data  pre-processing 

48  bad  channels  were  manually  identibed  by  visual  inspection.  After  discarding  these,  the 
data  set  contained  M  =  976  acquisitions.  Each  acquisition  contained  Nrec  =  976  time  traces. 
Each  time  trace  contained  L  —  2112  time  samples.  The  976  good  channels  were  indexed  from 
0  to  975.  The  corresponding  data  acquisitions  were  indexed  in  the  same  way.  A  Hann-window 
low-pass  filter  with  a  cutoff  frequency  of  4  MHz  was  applied  to  every  time  trace  in  both  the 
calibration  and  the  measurement  data.  This  data  filtering  was  implemented  to  mitigate  numerical 
errors  that  could  be  introduced  by  our  second-order  numerical  solver. 

C.  Estimation  of  excitation  pulse 

The  shape  of  the  excitation  pulse  was  estimated  as  the  time  trace  of  the  calibration  data  (after 
pre-processing)  received  by  the  488-th  receiver  at  the  0-th  data  acquisition.  Note  that  the  488-th 
receiver  was  approximated  located  on  the  axis  of  the  0-th  emitter,  thus  the  received  pulse  was 
minimally  affected  by  the  finite  aperture  size  effect  of  the  transducers.  Because  our  calibration 
data  and  measurement  data  were  acquired  using  different  electronic  amplifier  gains,  the  amplitude 
of  the  excitation  pulse  was  estimated  from  the  measurement  data.  More  specifically,  we  simulated 
the  0-th  data  acquisition  using  our  second-order  numerical  solver  and  compared  the  simulated 
time  trace  received  by  the  300-th  receiver  with  the  corresponding  measured  time  trace  (after  pre¬ 
processing).  The  ratio  between  the  maximum  values  of  these  two  traces  was  used  to  scale  the 
excitation  pulse  shape.  We  selected  the  300-th  receiver  because  it  resided  out  of  the  fan-region 
indicated  in  Fig.  1;  its  received  signals  were  unlikely  to  be  strongly  affected  by  the  presence  of 
the  object.  The  estimated  excitation  pulse  and  its  amplitude  spectrum  are  displayed  in  Fig.  16. 
Note  that  the  experimental  excitation  pulse  contained  higher  frequency  components  than  did  the 
computer-simulated  excitation  pulse  shown  in  Fig.  3. 

D.  Synthesis  of  combined  data 

As  discussed  in  Section  IV-C4,  signals  received  by  receivers  located  near  the  emitter  can  be 
unreliable  [23].  Our  experimental  data,  as  shown  in  Fig.  17-(a),  contained  noise-like  measure¬ 
ments  for  the  receivers  indexed  from  0  to  200,  and  from  955  to  975,  in  the  case  where  the  0-th 
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transducer  functioned  as  the  emitter.  Also,  our  point-like  transducer  assumption  introduces  larger 
model  mismatches  for  the  receivers  located  near  the  emitter.  As  shown  in  Figs.  17-(c)  and  -(d), 
even  though  the  simulated  time  trace  received  by  the  300-th  receiver  matches  accurately  with 
the  experimentally  measured  one,  the  simulated  time  trace  received  by  the  200-th  receiver  is 
substantially  different  compared  with  the  experimentally  measured  one.  In  order  to  minimize  the 
effects  of  model  mismatch,  we  replaced  these  unreliable  measurements  with  computer-simulated 
water  bath  data,  as  described  in  Section  IV-C.  We  designated  the  time  traces  received  by  the 
512  receivers  located  on  the  opposite  side  of  the  emitter  as  the  reliable  measurements  for  each 
data  acquisition.  The  0-th  data  acquisition  of  the  combined  data  is  displayed  in  Fig.  17-(b). 

E.  Estimation  of  initial  guess 

The  initial  guess  for  the  WISE  method  was  obtained  by  use  of  the  bent-ray  reconstruction 
method  described  in  Appendix  C.  We  first  filtered  each  time  trace  of  the  raw  data  by  a  band¬ 
pass  butterworth  filter  (0.5MHz  -  2.5MHz).  Subsequently,  we  extracted  the  TOF  by  use  of  the 
thresholding  method  with  a  thresholding  value  of  20%  of  the  peak  value  of  each  time  trace. 
The  bent-ray  reconstruction  algorithm  was  applied  for  image  reconstruction  with  a  measured 
background  sound  speed  1.513  mm//is.  The  resulting  image  is  shown  in  Fig.  18-(a)  and  has  a 
pixel  size  of  1  mm.  Finally,  the  image  was  smoothed  by  convolving  it  with  a  2D  Gaussian  kernel 
with  a  standard  deviation  of  2  mm. 

F.  Image  reconstruction 

We  applied  the  WISE  method  with  a  TV  penalty  to  the  combined  data  set.  The  second-order 
numerical  solver  was  employed  with  a  calculation  domain  of  dimensions  512.0  x  512.0  mm2. 
The  calculation  domain  was  sampled  on  a  2560  x  2560  Cartesian  grid  with  a  grid  spacing  of 
0.2  mm.  On  a  platform  consisting  of  dual  quad-core  CPUs  with  a  3.30  GHz  clock  speed,  64 
GB  RAM,  and  a  single  NVIDIA  Tesla  K20  GPU,  each  numerical  solver  run,  took  40  seconds 
to  calculate  the  pressure  data  for  2112  time  samples.  Knowing  the  size  of  the  phantom,  we  set 
the  reconstruction  region  to  be  within  a  circle  of  diameter  128  mm,  i.e.,  only  the  sound  speed 
values  of  pixels  within  the  circle  were  updated  during  the  iterative  image  reconstruction.  We 
swept  the  value  of  /3TV  over  a  wide  range  to  investigate  its  impact  on  the  reconstructed  images. 
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G.  Images  reconstructed  from  experimental  data 

As  shown  in  Fig.  18,  the  spatial  resolution  of  the  image  reconstructed  by  use  of  the  WISE 
method  with  a  TV  penalty  is  significantly  higher  than  that  reconstructed  by  use  of  the  bent-ray 
model-based  method.  In  particular,  the  structures  labeled  ‘A’  and  ‘B’  possess  clearly-defined 
boundaries.  This  observation  is  further  confirmed  by  the  profiles  of  the  two  images  shown  in 
Fig.  19.  In  addition,  the  structure  labeled  ‘C’  in  Fig.  18-(b)  is  almost  indistinguishable  in  the 
image  reconstructed  by  use  of  the  bent-ray  model-based  method  (see  Fig.  18-(a)).  The  improved 
spatial  resolution  is  expected  because  the  WISE  method  with  a  TV  penalty  takes  into  account  the 
high-order  diffraction,  which  is  ignored  by  the  bent-ray  method  [23].  Though  not  shown  here, 
for  the  bent-ray  method,  we  investigated  multiple  time-of-flight  pickers  [25],  and  systematically 
tuned  the  regularization  parameter.  As  such,  it  is  likely  that  Fig.  18-(a)  represents  a  nearly  optimal 
bent-ray  image  in  terms  of  the  resolution.  This  resolution  also  appears  to  be  similar  to  previous 
experimental  results  reported  in  the  literature  [26]. 

The  convergence  properties  of  the  WISE  method  with  a  TV  penalty  with  experimental  data 
were  consistent  with  those  observed  in  the  computer-simulation  studies.  Images  reconstructed  by 
use  of  10,  50,  and  300  algorithm  iterations  are  displayed  in  Fig.  20.  The  image  reconstructed  by 
use  of  10  iterations  contains  radial  streak  artifacts  that  are  similar  in  nature  to  those  observed  in 
the  computer- simulation  studies.  These  artifacts  were  mitigated  after  more  iterations.  The  image 
reconstructed  after  300  iterations  (Fig.  20-(d))  appears  to  be  similar  to  that  after  200  iterations 
(Fig.  18-(b)),  suggesting  that  the  WISE  method  with  a  TV  penalty  is  close  to  convergence  after 
about  200  iterations.  The  computational  time  for  completing  200  iterations  was  approximately 
14  hours.  The  estimated  computational  time  for  the  conventional  method  was  about  one  month, 
assuming  the  same  number  of  iterations  is  required  as  in  the  computer-simulation  studies  (i.e., 
40). 

Despite  the  nonlinearity  of  the  WISE  method  with  a  TV  penalty,  the  impact  of  the  TV 
smoothness  penalty  appears  to  be  similar  to  that  observed  in  other  imaging  applications  [51], 
[54]  (see  Fig.  21).  Though  not  shown  here,  the  impact  of  the  quadratic  penalty  is  also  similar. 
As  expected,  a  larger  value  of  /3  reduced  the  noise  level  at  the  expense  of  spatial  image 
resolution.  These  results  suggest  a  predictable  impact  of  the  smoothness  penalties  on  the  images 
reconstructed  using  the  WISE  method  with  a  TV  penalty. 
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VII.  Summary 

It  is  known  that  waveform  inversion-based  reconstruction  methods  can  produce  sound  speed 
images  that  possess  improved  spatial  resolution  properties  over  those  produced  by  ray-based 
methods.  However,  waveform  inversion  methods  are  computationally  demanding  and  have  not 
been  applied  widely  in  USCT  breast  imaging.  In  this  work,  based  on  the  time-domain  wave 
equation  and  motivated  by  recent  mathematical  results  in  the  geophysics  literature,  the  WISE 
method  was  developed  that  circumvents  the  large  computational  burden  of  conventional  wave¬ 
form  inversion  methods.  This  method  encodes  the  measurement  data  using  a  random  encoding 
vector  and  determines  an  estimate  of  the  sound  speed  distribution  by  solving  a  stochastic  opti¬ 
mization  problem  by  use  of  a  stochastic  gradient  descent  algorithm.  With  our  current  GPU-based 
implementation,  the  computation  time  was  reduced  from  weeks  to  hours.  The  WISE  method  was 
systematically  investigated  in  computer-simulation  and  experimental  studies  involving  a  breast 
phantom.  The  results  suggest  that  the  method  holds  value  for  USCT  breast  imaging  applications 
in  a  practical  setting. 

Many  opportunities  remain  to  further  improve  the  performance  of  the  WISE  method.  As  shown 
in  Fig.  18,  images  reconstructed  by  use  of  the  WISE  method  can  contain  certain  artifacts  that 
are  not  present  in  the  image  reconstructed  by  use  of  the  bent-ray  method.  An  example  of  such 
an  artifact  is  the  dark  horizontal  streak  below  the  structure  C.  Because  of  the  nonlinearilty  of  the 
image  reconstruction  problem,  it  is  challenging  to  determine  whether  these  artifacts  are  caused 
by  imaging  model  errors  or  by  the  optimization  algorithm,  which  might  have  arrived  at  a  local 
minimum  of  the  cost  function.  A  more  accurate  imaging  model  can  be  developed  to  account  for 
the  out-of-plane  scattering,  the  transducer  finite  aperture  size  effect,  and  the  medium  acoustic 
absorption,  etc.  Also,  the  stochastic  gradient  descent  algorithm  is  one  of  the  most  basic  stochastic 
optimization  algorithms.  Numerous  emerging  optimization  algorithms  can  be  employed  [33],  [41] 
to  improve  the  convergence  rate,  as  well  as  to  avoid  local  minima. 

There  remains  a  need  to  conduct  additional  investigations  of  the  numerical  properties  of  the 
WISE  method.  Currently,  a  systematic  comparsion  of  the  statistical  properties  of  the  WISE  and 
the  sequential  waveform  inversion  methods  is  prohibted  by  the  3-day  computational  time  for 
the  latter  method.  This  comparsion  will  be  interesting  when  a  more  efficient  numerical  solver  is 
available.  Given  the  fact  that  waveform  inversion  is  nonlinear  and  sensitive  to  its  initial  guess, 
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it  becomes  important  to  investigate  how  to  obtain  an  accurate  and  computationally  efficient 
initial  guess.  We  also  observed  that  the  performance  of  the  WISE  method  is  sensitive  to  how 
strong  the  medium  heterogeneities  are  and  the  profile  of  the  excitation  pulse.  Investigation  of  the 
impact  of  the  excitation  pulse  on  the  image  reconstruction  may  help  optimize  hardware  design. 
In  addition,  exploring  the  statistics  of  the  reconstructed  images  will  allow  a  quantitative  image 
quality  assessment  for  cancer  diagnosis. 


Appendix  A 

Continuous-to-Discrete  USCT  Imaging  Model 


In  practice,  each  data  function  gm(r,t)  is  spatially  and  temporally  sampled  to  form  a  data 
vector  gm  e  M;V"  ,  L,  where  NTec  and  L  denote  the  number  of  receivers  and  the  number  of  time 
samples,  respectively.  We  will  assume  that  Nrec  and  L  do  not  vary  with  excitation  pulse.  Let 
[gm]„reci+|  denotes  the  (nrecL  +  /)-th  element  of  gm.  When  the  receivers  are  point-like,  gm  is 
defined  as 

[Sm\niecL+i  =  9m{r{m ,  nrec),  l A1),  (29) 


where  the  indices  nrec  and  /  specify  the  receiver  location  and  temporal  sample,  respectively,  and 
A*  is  the  temporal  sampling  interval.  The  vector  r(m,  nrec)  e  1  lrn  denotes  the  location  of  the 
nrec-th  receiver  at  the  m-th  data  acquisition. 

A  C-D  imaging  model  for  USCT  describes  the  mapping  of  c(r)  to  the  data  vector  gm  and 
can  be  expressed  as 


[Sm]niecL+(  AA.mH,  Sm{v,t') 


r=r(m,nrec),i=/At 


for 


nrec=0,l,-,iVrec-l 
1=0,1,—  ,L- 1 


(30) 


Note  that  the  acousto-electrical  impulse  response  of  the  receivers  can  be  incorporated  into  the 
C-D  imaging  model  by  temporally  convolving  sm(r,t)  in  Eqn.  (1)  with  the  receivers’  acousto¬ 
electrical  impulse  response  if  we  assume  all  receiving  transducers  share  an  identical  acousto¬ 
electrical  impulse  response. 


Appendix  B 

Frechet  derivative  of  data  fidelity  term 
Consider  the  integrated  squared-error  data  misfit  function,  [22],  [23] 

M~1  p  pT 

•7rCC(c)  =  x  ^2  /  dr  dt\gm(r,t)  -  gm(r,t)]2,  (31) 

m= 0 


September  26,  2014 


DRAFT 


25 


where  gm(r,t)  and  gm(r,t)  denote  the  measured  data  function  and  the  predicted  data  function 
computed  by  use  of  Eqn.  (3)  with  the  current  estimate  of  c(r). 

Both  the  sequential  and  WISE  reconstruction  method  described  in  Section  III  require  knowl¬ 
edge  of  the  Frechet  derivatives  of  T(  (  {c)  and  lZcc(c)  with  respect  to  c,  denoted  by  Vc  J700 
and  Vc7£cc,  respectively.  The  calculation  of  Vc7 Zcc  can  be  readily  accomplished  for  quadratic 
smoothness  penalties  [51],  [55].  For  the  integrated  squared  error  data  misfit  function  given  in 
Eqn.  (31),  Vc J700  can  be  computed  via  an  adjoint  state  method  as  [28],  [56],  [57] 

1  M_1  fT  d2 

VcJrCC  =  rTH^/  dtVm{r,T-t)—pm(r,t),  (32) 

where  qm(r,t)  G  L2(E3  x  [0,  oo))  is  the  solution  to  the  adjoint  wave  equation.  The  adjoint  wave 
equation  is  defined  as 


V2<?m(r,f)  - 


i  d2 


■qm(r,t)  =  —Tm(r,  t), 


(33) 


c2(r) d2t 

where  rm(r,  t)  =  gm(r,T  —  t)  —  gm(r,T  —  t ).  The  adjoint  wave  equation  is  nearly  identical  in 
form  to  the  wave  equation  in  Eqn.  (1)  except  for  the  different  source  term  on  the  right-hand 
side,  suggesting  the  same  numerical  approach  can  be  employed  to  solve  both  equations.  Since 
one  needs  to  solve  Eqns.  (1)  and  (33)  M  times  in  order  to  calculate  VcJr<  c,  it  is  generally  true 
that  the  sequential  waveform  inversion  is  computationally  demanding  even  for  a  2D  geometry 


[58]. 


Appendix  C 

Bent-ray  model-based  sound  speed  reconstruction 

We  developed  an  iterative  image  reconstruction  algorithm  based  on  a  bent-ray  imaging  model. 
The  bent-ray  imaging  model  assumes  that  an  acoustic  pulse  travels  along  a  ray  path  that  connects 
the  emitter  and  the  receiver  and  accounts  for  the  refraction  of  rays,  also  known  as  ray-bending, 
through  an  acoustically  inhomogeneous  medium.  For  each  pair  of  receiver  and  emitter,  the  travel 
time,  as  well  as  the  ray  path,  is  determined  by  the  medium’s  sound  speed  distribution.  Given 
the  travel  times  for  a  collection  of  emitter- and-receiver  pairs  distributed  around  the  object,  the 
medium  sound  speed  distribution  can  be  iteratively  reconstructed.  This  bent-ray  model-based 
sound  speed  reconstruction  (BRSR)  method  has  been  employed  in  the  USCT  literature  [26], 
[59],  [60]. 
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In  order  to  perform  the  BRSR,  we  extracted  a  TOF  data  vector  from  the  measured  pressure 
data.  Denoting  the  TOF  data  vector  by  T  e  MMArrec,  each  element  of  T  represented  the  TOF 
from  each  emitter-and-receiver  pair.  The  extraction  of  the  TOF  was  conducted  in  two  steps.  First, 
we  estimated  the  difference  between  the  TOF  when  the  object  was  present  and  the  TOF  when 
the  object  was  absent  by  use  of  a  thresholding  method  [25],  [61].  In  particular,  20%  of  the  peak 
value  of  each  time  trace  was  employed  as  the  thresholding  value.  Second,  a  TOF  offset  was 
added  to  the  estimated  difference  TOF  for  each  emitter-and-receiver  pair  to  obtain  the  absolute 
TOF,  where  the  TOF  offset  was  calculated  according  to  the  scanning  geometry  and  the  known 
background  SOS. 

Having  the  TOF  vector  T,  we  reconstructed  the  sound  speed  by  solving  the  following  opti¬ 
mization  problem: 

s  =  argmin  ||  T  —  Kss  ||2  +/37l(s),  (34) 

S 

where  s  denotes  the  slowness  (the  reciprocal  of  the  SOS)  vector,  and  Ks  denotes  the  system 
matrix  that  maps  the  slowness  distribution  to  the  TOF  data.  The  superscript  ‘s’  indicates  the 
dependence  of  Ks  on  the  slowness  map.  At  each  iteration,  using  the  current  estimate  of  the  SOS, 
a  ray-tracing  method  [62]  was  employed  to  construct  the  system  matrix  Ks.  Explicitly  storing 
the  system  matrix  in  the  sparse  representation,  we  utilized  the  limited  BFGS  method  [63]  to 
solve  the  optimization  problem  given  in  Eqn.  (34).  The  estimated  slowness  was  then  converted 
to  the  sound  speed  by  taking  the  reciprocal  of  s  element- wisely.  We  refer  the  readers  to  [26], 
[59]— [61],  [64]  for  more  details  about  the  BRSR  method. 
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Tables 

TABLE  I:  Parameters  of  the  numerical  breast  phantom 


Structure 

index 

Tissue  type 

Sound  speed 

[mm-/xs_1] 

Slope  of  attenuation 
[dB  ■  (MHz)  ~y  -cm~ 1  ] 

0 

Adipose 

1.47 

0.60 

1 

Parenchyma 

1.51 

0.75 

2 

Benign  tumor 

1.47 

0.60 

3 

Benign  tumor 

1.47 

0.60 

4 

Cyst 

1.53 

0.00217 

5 

Malignant  tumor 

1.565 

0.57 

6 

Malignant  tumor 

1.565 

0.57 

7 

Malignant  tumor 

1.57 

0.57 
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TABLE  II:  Parameters  of  the  experimental  breast  phantom 


Material 

Sound  speed 

[irmi'/is-1] 

Attenuation  coefficient 

at  2.5  MHz  [dB/cm] 

Fat 

1.467 

0.48 

Parenchymal  tissue 

1.552 

0.89 

Cancer 

1.563 

1.20 

Fibroadenoma 

1.552 

0.52 

Gelatin  cyst 

1.585 

0.16 
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Figures 


Fig.  1:  Schematic  of  a  USCT  system  with  a  circular  transducer  array  whose  elements  are 
indexed  from  0  to  255.  It  shows  the  first  data  acquisition,  where  element-0  (in  red)  is  emitting 
an  acoustic  pulse,  while  all  256  elements  are  receiving  signals.  The  region-of-interest  (ROI)  is 
shaded  in  gray,  and  the  dashed  square  box  represents  the  physical  dimensions  (128  x  128  mm2) 
of  all  reconstructed  images. 
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Fig.  2:  (a)  Sound  speed  map  [mm-/is  []  and  (b)  acoustic  absorption  coefficient  map 

[dB-(MHz)_y-cm~1]  of  the  numerical  breast  phantom. 
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Fig.  3:  (a)  Normalized  temporal  profile  and  (b)  amplitude  spectrum  of  the  excitation  pulse 
employed  in  the  computer-simulation  studies.  The  dashed  line  in  (b)  marks  the  center  frequency 
of  excitation  pulse  at  0.82  MHz. 
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Fig.  4:  Computer-simulated  (a)  noise-free  and  (b)  noisy  data  vectors  at  the  0-th  data  acquisition, 
(c)  Profiles  of  the  pressure  received  by  the  128-th  transdcuer.  The  grayscale  window  for  (a)  and 
(b)  is  [-45,  0]  dB. 
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(a)  (b) 


(c)  (d) 

Fig.  5:  Images  reconstructed  by  use  of  (a)  the  WISE  method  after  the  199-th  iteration  (1,  018  runs 
of  the  numerical  solver),  (b)  the  sequential  waveform  inversion  algorithm  after  the  43-rd  iteration 
(57,  088  runs  of  the  numerical  solver),  (c)  the  bent-ray  model-based  sound  speed  reconstruction 
method,  and  (d)  the  sequential  waveform  inversion  algorithm  after  the  1-st  iteration  (1,024 
runs  of  the  numerical  solver)  from  the  noise-free  non-attenuated  data.  The  grayscale  window  is 
[1.46, 1.58]  mm//is. 
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Fig.  6:  Profiles  at  y  =  6.5  mm  of  the  images  reconstructed  by  use  of  the  bent-ray  TOF  image 
reconstruction  method  and  the  WISE  method  from  the  noise-free  non-attenuated  data. 
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(c)  (d) 

Fig.  7:  Images  reconstructed  by  use  of  the  WISE  method  after  (a)  the  20-th,  (b)  the  50-th,  (c)  the 
100-th,  and  (d)  the  250-th  iteration  from  the  noise-free,  non-attenuated  data  set.  The  grayscale 
window  is  [1.46, 1.58]  mm//is. 
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Fig.  8:  Profiles  of  the  images  reconstructed  by  use  of  the  WISE  method  from  the  noise-free 
non-attenuated  data  after  different  numbers  of  iterations. 
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Fig.  9:  Plots  of  the  absolute  Euclidean  distances  versus  (a)  the  number  of  iterations  and  (b) 
the  number  of  numerical  solver  runs,  (c)  Plots  of  the  cost  function  value  versus  the  number  of 
iterations. 
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Fig.  10:  Images  reconstructed  from  non-attenuated  data  contaminated  with  Gaussian  random 
noise.  Images  (a-c)  were  reconstructed  by  use  of  the  WISE  method  with  a  quadratic  penalty 
with  ft®  =  1.0  x  10'3,  1.0  x  10-2,  and  1.0  x  10_1,  respectively.  Image  (d)  was  reconstructed  by 
use  of  the  WISE  method  with  a  TV  penalty  with  /3TV  =  5.0  x  10~4.  The  insert  in  the  up  right 
corner  of  each  image  is  the  zoomed-in  image  of  the  dashed  black  box.  The  grayscale  window 
is  [1.46, 1.58]  mm//is. 
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Fig.  11:  (a)  Computer- simulated  noise-free  attenuated  pressure  of  the  0-th  data  acquisition. 

(b)  The  difference  between  the  attenuated  pressure  data  and  the  non-attenuated  pressure  data. 

(c)  The  temporal  profiles  and  (d)  the  amplitude  spectra  of  the  pressure  received  by  the  128-th 
transducer.  The  grayscale  window  for  (a)  and  (b)  is  [—45,0]  dB. 
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(a)  (b) 


Fig.  12:  (a)  Image  reconstructed  by  use  of  the  WISE  method  from  the  noise-free  attenuated  data, 
(b)  Image  reconstructed  by  use  of  the  WISE  method  with  a  TV  penalty  with  /3TV  =  5.0  x  HE4, 
from  the  noisy  attenuated  data.  The  grayscale  window  is  [1.46, 1.58]  mm//is.  (c)  Profiles  at 
y  —  6.5  mm  of  the  reconstructed  images. 
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Fig.  13:  (a)  Computer-simulated  noise-free  non-attenuated  pressure  data  when  the  object  is 

absent,  (b)  The  difference  between  the  pressure  data  when  object  is  present  and  the  pressure 
data  when  the  object  is  absent,  (c)  Profiles  of  the  pressure  received  by  the  40-th  transducer.  The 
grayscale  window  for  (a)  and  (b)  is  [—45,  0]  dB. 
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(a)  (b) 


Fig.  14:  Images  reconstructed  by  use  of  the  WISE  method  from  noise-free  combined  data  that 
are  completed  (a)  with  computer-simulated  pressure  corresponding  to  a  homogeneous  medium 
and  (b)  with  zeros.  The  grayscale  window  is  [1.46, 1.58]  mm//is.  (c)  Profiles  at  y  =  6.5  mm  of 
the  images  reconstructed  by  use  of  the  WISE  method  from  the  two  combined  data  sets. 
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Fibroadenoma 


Fig.  15:  Schematic  of  the  breast  phantom  employed  in  the  experimental  study. 


September  26,  2014 


DRAFT 


48 


(a)  (b) 

Fig.  16:  (a)  Normalized  temporal  profile  and  (b)  amplitude  spectrum  of  the  excitation  pulse 
employed  in  the  experimental  studies.  The  dashed  line  in  (b)  marks  the  center  frequency  of 
excitation  pulse  at  2.09  MHz. 
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Fig.  17:  Zeroth  acquisition  of  (a)  the  experimentally-measured  raw  data  and  (b)  the  combined 
data,  respectively,  and  time  traces  at  the  0-th  acquisition  received  by  (c)  the  300-th  receiver,  and 
(d)  the  200-th  receiver,  respectively.  The  grayscale  window  for  (a)  and  (b)  is  [—45,  0]  dB. 
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(a)  (b) 

Fig.  18:  Images  reconstructed  from  the  experimentally  measured  phantom  data  by  use  of  (a)  the 
bent-ray  model-based  sound  speed  reconstruction  method  and  (b)  the  WISE  method  with  a  TV 
penalty  with  (/3TV  =  1.0  x  102)  after  the  200-th  iteration.  The  grayscale  window  is  [1.49, 1.57] 
mm//is. 
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Fig.  19:  Profiles  at  (a)  x  =  —24.0  mm  and  (b)  x  =  10.0  mm  of  the  reconstructed  images  by 
use  of  the  bent-ray  model-based  sound  speed  reconstruction  method  (light  solid)  and  the  WISE 
method  with  a  TV  penalty  with  /3TV  =  1.0  x  102  (dark  dashed)  from  experimentally  measured 
data. 
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(c)  (d) 

Fig.  20:  (a)  The  initial  guess  of  the  sound  speed  map  and  the  images  reconstructed  by  use  of  the 
WISE  method  with  a  TV  penalty  with  (/3TV  =  f.O  x  102)  after  (b)  the  10-th,  (b)  the  50-th  and 
(d)  the  300-th  iteration,  from  the  experimentally  measured  phantom  data.  The  grayscale  window 
is  [1.49, 1.57]  mm//is. 
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(a)  (b) 

Fig.  21:  Images  reconstructed  by  use  of  the  WISE  method  with  a  TV  penalty  with  (a)  /3TV  = 
5.0  x  101,  and  (b)  /3TV  =  5.0  x  102,  from  the  experimentally  measured  phantom  data.  The 
grayscale  window  is  [1.49, 1.57]  mm//is. 
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